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Abstract 


In  this  thesis  we  present  results  from  inversion  of  data  using  dense  arrays  of 
collocated  seismic  and  magnetotelluric  stations  located  in  the  Cascadia 
subduction  zone  region  of  central  Washington.  In  the  migrated  seismic  section, 
we  clearly  image  the  top  of  the  slab  and  oceanic  Moho,  as  well  as  a  velocity 
increase  corresponding  to  the  eclogitization  of  the  hydrated  upper  crust.  A 
deeper  velocity  increase  is  interpreted  as  the  eclogitization  of  metastable 
gabbros,  assisted  by  fluids  released  from  the  dehydration  of  upper  mantle 
chlorite.  A  low  velocity  feature  interpreted  as  a  fluid/melt  phase  is  present  above 
this  transition.  The  serpentinized  wedge  and  continental  Moho  are  also  imaged. 
The  magnetotelluric  image  further  constrains  the  fluid/melt  features,  showing  a 
rising  conductive  feature  that  forms  a  column  up  to  a  conductor  indicative  of  a 
magma  chamber  feeding  Mt.  Rainier.  This  feature  also  explains  the  disruption  of 
the  continental  Moho  found  in  the  migrated  image.  Exploration  of  the 
assumption  of  smoothness  implicit  in  the  standard  MT  inversion  provides  tools 
that  enable  us  to  generate  a  more  accurate  MT  model.  This  final  MT  model 
clearly  demonstrates  the  link  between  slab  derived  fluids/melting  and  the  Mt. 
Rainier  magma  chamber. 
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1.  Introduction 


Many  of  the  natural  hazards  associated  with  subduction  zones  are  intimately 
tied  to  the  release  of  fluids  bound  within  the  subducting  slab.  Intraslab 
earthquakes  have  been  linked  with  a  process  in  which  fluids  released  during  the 
transformation  of  hydrated  metabasalts  to  eclogite  in  the  upper  crust  increases 
pore  pressure,  thereby  reactivating  pre-existing  faults  (Kirby  et  al.,  1996,  Hacker 
et  al.,  2003b,  Preston  et  al.,  2003).  Slow  slip  and  low  frequency  tremor,  which 
have  been  tentatively  linked  with  the  periodicity  of  great  megathrust 
earthquakes  (Mazzotti  and  Adams,  2004),  can  be  explained  by  episodic  buildup 
and  release  of  fluid  pore  pressure  across  plate  boundaries  down-dip  of  the 
locked  zone  (Abers  et  al.,  2009,  Audet  et  al.,  2010).  Most  arc  magmas  are  also 
generally  believed  to  be  generated  by  flux  melting  that  occurs  when  dehydration 
reactions  in  the  oceanic  crust  and  uppermost  mantle  release  fluids  into  the 
continental  asthenosphere  (Ulmer  and  Trommsdorff,  1995,  Kirby  et  al.,  1996). 
Finally,  density  increases  associated  with  dehydration  have  been  linked  with 
changes  in  slab  dip  and  increased  stresses  within  the  slab  (Klemd  et  al.,  2011). 
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While  recent  investigations  have  dramatically  improved  our  understanding  of 
the  distribution  and  cycling  of  fluids  in  subduction  systems,  many  questions 
remain  open.  Advances  in  the  understanding  of  the  phase  transformations  that 
occur  within  the  subducting  lithosphere  (Poli  and  Schmidt,  2002,  Hacker  et  al., 
2003a  and  references  within,  Hacker,  2008)  have  given  us  the  tools  to  better 
interpret  geophysical  evidence.  Improvements  in  the  thermal  modeling  of 
subduction  zones  such  as  the  use  of  olivine  rheology  in  the  mantle  wedge  (van 
Keken  et  al.,  2002)  and  decoupling  at  the  cold  nose  of  the  wedge  (Wada  and 
Wang,  2009)  have  aided  in  the  development  of  a  consistent  set  of  thermal  models 
for  subduction  systems  worldwide  (Syracuse  et  al.,  2010).  Application  of  the 
phase  transformations  along  the  P-T  regimes  predicted  by  the  thermal  models 
has  made  it  possible  to  predict  depth  dependent  fluid  flux  for  the  global  suite  of 
subduction  zones  (van  Keken  et  al.,  2011). 

Cascadia  is  an  attractive  setting  for  studying  the  fluid  distribution  and 
dynamics  in  subduction  zones  for  a  variety  of  reasons.  It  has  been  thoroughly 
investigated  using  a  variety  of  geological  and  geophysical  methods  over  a 
number  of  decades,  so  much  is  known  about  both  the  surface  geology  and  the 
structure  of  the  subsurface.  As  the  slab  is  young  and  hot,  it  is  expected  to  release 
a  significant  fraction  of  its  fluids  in  the  forearc,  suggesting  that  the  evidence  of 
fluids  should  be  more  pronounced  here  than  they  would  be  in  a  cooler  setting. 
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Finally,  the  connection  between  these  fluid  processes  and  the  associated  tectonic 
and  volcanic  hazards  that  potentially  threaten  a  number  of  densely  populated 
areas  underline  the  importance  of  understanding  the  subduction  fluid  cycle  in 
Cascadia. 

For  all  of  these  reasons,  the  Cascadia  region  has  been  the  subject  of  extensive 
geophysical  work.  While  much  of  the  geophysical  exploration  of  the  region  has 
involved  seismology,  the  magnetotelluric  method  has  also  been  employed  to  good 
effect.  Magnetotellurics  is  attractive  because  while  it  is  sensitive  to  a  completely 
different  set  of  rock  properties  than  are  seismic  methods  (Hellfrich,  2003,  Jones  et  al., 
2009),  velocity  and  conductivity  are  often  correlated  (Stanley  et  al.,  1990,  Jones  et  al., 
2009),  either  directly  or  simply  in  the  fact  that  at  distinct  geological  boundaries  both 
parameters  are  likely  to  change  (Gallardo,  2007,  Moorkamp  et  al.,  2007).  It  is  therefore 
tempting  to  utilize  both  types  of  data  in  our  interpretations,  particularly  where  the  data 
profiles  are  collocated  (Gallardo  and  Meju,  2007). 

While  much  work  in  the  past  has  simply  been  a  matter  of  questioning  whether  an 
interpretation  of  one  type  of  data  is  consistent  with  a  previous  interpretation  of  a 
different  type,  there  have  been  several  recent  efforts  to  combine  two  (or  more)  types  of 
data  on  a  more  fundamental  level.  The  rationale  behind  joint  inversion  of  multiple  data 
sets  is  to  improve  resolution,  diminish  the  influence  of  noise,  and  reduce  the  inherent 
non-uniqueness  by  providing  additional  constraints  (Gallardo,  2007,  Moorkamp  et  al.. 
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2007).  While  the  values  of  velocity  and  resistivity  have  been  compiled  for  a  wide  range 
of  geological  and  physical  conditions  (Hacker,  2003a,  Jones  et  al.,  2009),  there  exists  no 
fundamental  relationship  between  the  two  physical  properties  (Bedrosian  et  al.,  2007). 
Still,  as  we  stated  above,  the  two  properties  do  tend  to  be  reasonably  well  correlated  in 
the  subsurface,  and  several  different  approaches  to  joint  inversion  of  seismic  and 
electromagnetic  data  have  been  attempted.  Some  of  the  techniques  utilized  to  date 
include  stochastic  optimization  methods  such  as  genetic  algorithms  (Moorkamp  et  al., 
2007),  physical  property  gradient  methods  (often  combined  with  structural  constraints) 
such  as  the  cross-gradient  (Gallardo,  2007,  Hu  et  al.,  2008,  Gallardo  and  Meju,  2007), 
and  statistical  correlation  methods  such  as  cluster  analysis  (Bedrosian  et  al.,  2007, 
Munoz  et  al.,  2010). 

None  of  the  joint  inversion  methods  mentioned  above  address  this  particular 
combination  of  ideas  in  a  completely  satisfactory  way.  They  all  rely  in  some  way  on  the 
natural  relationship  between  the  distribution  of  the  resistivity  and  velocity  values.  The 
distribution  of  resistivity  values  near  sharp  boundaries  in  the  MT  inversion  is  smoothed 
by  design  (which  directly  impacts  the  gradient  as  well),  creating  an  artificial  difference 
that  the  joint  inversions  would  attempt  to  resolve  as  a  natural  phenomenon. 

A  different  way  to  look  at  the  use  of  multiple  data  sets  is  to  attempt  to  use  the 
strengths  of  one  technique  to  mitigate  weaknesses  inherent  in  another.  For  example, 
the  magnetotelluric  technique  typically  employs  a  smoothness  constraint  to  control  the 
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resolution  of  structure  to  that  required  by  the  data.  Of  course,  sharp  boundaries  exist  in 
the  subsurface,  and  these  are  often  not  well  resolved  by  the  method.  This  suggests  that 
we  might  attempt  to  combine  MT  with  a  method  that  is  well  suited  to  illuminating 
sharp  boundaries.  Teleseismic  migration  is  one  particularly  appealing  candidate, 
involving  backpropagation  of  a  scattered  wavefield  to  the  scattering  points  in  the 
subsurface  (Rondenay,  2009),  thereby  highlighting  the  very  boundaries  we  hope  to 
resolve. 

This  is  the  approach  we  take  in  the  current  study,  analyzing  data  collected  at  dense 
arrays  of  roughly  collocated  passive  seismic  and  magnetotelluric  stations.  We  analyze 
the  two  data  sets  independently,  using  Generalized  Radon  Transform  (GRT)  migration 
for  the  seismic  data,  and  Non-Linear  Conjugate  Gradient  (NLCG)  inversion  for  the 
magnetotelluric  data.  We  then  perform  the  magnetotelluric  inversion  again,  this  time 
incorporating  the  a  priori  information  obtained  during  the  GRT  migration  into  the 
inversion. 

While  there  are  valid  reasons  for  searching  for  the  smoothest  resistivity  model  in  the 
absence  of  a  priori  information,  the  presence  of  such  information,  particularly  with 
respect  to  the  location  of  sharp  boundaries,  warrants  another  approach.  While 
inversion  methods  allow  us  to  create  models  that  fit  the  data  as  closely  as  desired,  there 
is  a  tradeoff  between  resolving  structure  and  generating  false  structure  from  noise  that 
tips  decidedly  in  favor  of  the  latter  as  the  misfit  becomes  small.  For  any  reasonable 
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selection  of  misfit,  there  are  an  infinite  number  of  solutions  that  fit  the  data  equally  well 
(Parker,  1994,  Zhdanov,  1993).  The  selection  of  a  specific  model  is  necessarily  a 
subjective  process  involving  the  incorporation  of  a  priori  information  and/or 
assumptions  about  the  nature  of  the  distribution  of  the  measured  parameter  in  the 
subsurface. 

The  conventional  way  to  think  about  such  problems  is  as  a  minimization  of  the 
Tikhonov  parametric  functional  (Tikhonov  and  Arsenin,  1977), 

T(m)=cj)(m)+As(m)  (1) 

where  m  is  a  matrix  representing  a  distribution  of  geophysical  parameters  within  some 
space,  cj)(m)  is  the  misfit  functional  described  by 

<j)(m)=  I  I  Fm-d  I  1 2  (2) 

where  d  is  the  observed  data  and  Fm  the  data  predicted  by  the  model,  and  s(m)  is  a 
functional  that  represents  assumptions  or  a  priori  information  that  the  worker  has 
chosen  to  include,  with  its  relative  importance  determined  by  the  regularization 
parameter  A  (Tikhonov  and  Arsenin,  1977,  Parker,  1994,  Mehanee  and  Zhdanov,  2002). 
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Equations  (1)  and  (2)  offer  at  least  three  fundamental  ways  in  which  we  might 
express  a  preference  for  one  model  over  another.  The  first  of  these  is  simply  to  fix  the 
value  of  the  model  parameter  over  some  spatial  distribution.  For  example,  in  coastal 
magnetotelluric  surveys,  it  is  standard  practice  to  include  in  the  model  a  fixed 
distribution  of  highly  conductive  cells  to  represent  the  ocean,  thereby  excluding  all 
models  that  don't  include  a  conductor  where  we  know  one  to  exist. 

A  second  method  is  to  include  an  additional  parameter  within  the  misfit  functional, 
jointly  inverting  for  two  nominally  independent  data  sets.  The  inclusion  of  a  second  set 
of  data  can  improve  resolution  and  reduce  the  impact  of  noise,  thereby  constraining  the 
range  of  possible  models  (Gallardo,  2007).  The  set  of  models  that  explains  two  (or 
more)  datasets  to  a  given  misfit  is  both  smaller  and  more  compelling  than  the  set  that 
explains  one  but  not  necessarily  the  other. 

Finally,  we  can  utilize  the  stabilizing  functional  s(m)  to  express  preference  for  a 
given  subset  of  models.  While  several  stabilizing  functions  have  been  successfully  used 
(least  squares,  minimum  norm  of  the  gradient,  distance  from  an  a  priori  model,  etc.)  in 
various  geophysical  applications  (Smith  et  al.,  1999,  Mehanee  and  Zhdanov,  2002,  de 
Groot-Hedlin  and  Constable,  2004),  many  magnetotelluric  algorithms  employ  the 
minimum  norm  of  the  Faplacian,  which  will  seek  the  smoothest  model  possible  that  fits 
the  data  (for  a  given  misfit).  The  argument  is  that  by  seeking  the  smoothest  model 
(again,  for  a  sufficiently  large  misfit),  only  structure  explicitly  required  by  the  data  will 


12 


be  resolved  (de  Groot-Hedlin  and  Constable,  1990).  In  nature,  sharp  geoelectric 
boundaries  do  of  course  exist,  and  smooth  algorithms  may  not  sufficiently  resolve  these 
boundaries  (Smith  et  al.,  1999,  de  Groot-Hedlin  and  Constable,  2004,  Mehanee  and 
Zhdanov,  2002). 

In  order  to  overcome  these  difficulties,  several  algorithms  have  been  developed  that 
allow  for  sharp  boundaries  by  making  different  global  assumptions  about  the 
distribution  of  the  modeled  parameter  in  the  subsurface,  for  example  by  favoring 
models  made  up  of  layers  (Smith  et  al.,  1999,  de  Groot-Hedlin  and  Constable,  2004)  or 
containing  primarily  blocky  structure  (Mehanee  and  Zhdanov,  2002).  Of  course,  each  of 
these  global  algorithms  present  trade-offs,  and  are  effective  largely  to  the  extent  that  the 
underlying  assumptions  are  accurate  throughout  the  modeled  space.  While  a  priori 
knowledge  regarding  the  location  and  nature  of  sharp  discontinuities  can  be  an 
extremely  useful  implement  in  resolving  the  associated  structure  (de  Groot-Hedlin  and 
Constable,  1990),  its  utility  in  resolving  a  sufficiently  complicated  subsurface  is 
dependent  on  the  extent  to  which  it  can  be  effectively  applied  locally. 

Inversion  packages  such  as  WinGlink  and  Occam  allow  for  the  elimination  of  the 
penalty  for  roughness  along  presumed  discontinuities  (de  Groot-Hedlin  and  Constable, 
1990,  Rodi  and  Mackie,  2001),  and  while  this  tool  has  been  used  to  good  effect  (Matsuno 
et  al.,  2010),  it  is  but  a  tentative  first  step  in  exploring  the  possibilities  for  spatially- 
dependent  stabilizing  functionals.  While  de  Groot-Hedlin  and  Constable  rightly  point 
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out  that  one  should  use  such  tools  cautiously  as  "incorrect  placement  of  boundaries 
may  produce  misleading  results",  we  must  also  recognize  that  their  use  is  simply 
trading  one  set  of  assumptions  for  another  and  we  should  be  cautious  in  any  event. 

Chapter  2  presents  high  resolution  images  from  2-D  GRT  migration  of  data 
collected  from  a  dense  array  of  41  broadband  seismometers  located  in  a  roughly  east  to 
west  transect  across  central  Washington  State.  We  see  a  dipping  low  velocity  layer  that 
we  interpret  to  be  the  descending  oceanic  crust,  underlain  by  a  positive  velocity 
transition  that  we  take  to  be  the  oceanic  Moho.  A  fading  velocity  contrast  suggests  the 
transformation  of  hydrated  basalt  in  the  upper  crust  to  eclogite,  but  also  shows  the 
retention  of  fluids  to  depths  of  80  km  or  more.  A  low  velocity  feature  above  the  slab  is 
consistent  with  a  fluid/melt  phase  released  at  depths  of  80-100  km. 

Previous  GRT  work  conducted  across  Vancouver  Island/British  Columbia  to  the 
north  and  Oregon  to  the  south  allow  us  to  make  comparisons  and  investigate  regional 
implications.  Earthquake  distribution  combined  with  our  results  suggest  that  the  fluid 
release  at  80  km  depth  in  central  Washington  is  due  to  the  dehydration  of  chlorite 
harzburgite  in  the  uppermost  mantle  and  subsequent  transformation  of  lower  crustal 
metastable  gabbros  to  eclogite.  The  extent  to  which  this  occurs  in  central  Washington 
appears  to  be  unique  for  Cascadia. 
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Chapter  3  follows  up  with  high  resolution  images  generated  from  a  non-linear 
conjugate  gradients  inversion  of  magnetotelluric  data  collected  from  60  broadband  and 
20  long  period  stations  along  a  line  collocated  with  the  seismic  array.  The  down-going 
slab  appears  as  a  large  resistive  feature,  though  the  surface  is  not  very  well  defined. 
Likewise,  there  appear  to  be  two  rising  conductive  features,  one  rising  from  a  depth  of 
roughly  40  km  that  appears  to  be  connected  with  the  hydrated  basalt  to  eclogite  fluid 
release,  and  a  second  one  starting  deeper  and  to  the  east  that  rises  upwards  and 
seawards,  collecting  just  below  and  to  the  west  of  Mt.  Rainier. 

The  first  part  of  Chapter  4  investigates  the  consequences  of  the  standard  inclusion  of 
a  smoothness  parameter  in  the  Tikhonov  equation  for  magnetotelluric  inversions. 
Through  a  series  of  models,  we  explore  ways  to  better  understand  the  impact  of  the 
smoothness  assumption  in  various  settings,  and  also  methods  to  mitigate  its  impact.  In 
the  second  part  of  Chapter  4  we  apply  those  lessons  to  the  Cascadia  Array  for 
Earthscope  Magnetotelluric  (CAFE  MT)  data  set,  incorporating  the  results  of  the  seismic 
migration  into  the  MT  inversion  sequence.  We  conclude  by  showing  that  the 
augmented  model  is  more  robust  than  the  standard  half-space  models,  and  therefore 
provides  stronger  insight  into  the  fluid  processes  of  the  Cascadia  subduction  zone  in 
Washington. 
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Abstract 

In  this  paper  we  present  high  resolution  seismic  images  from  2-D  Generalized 
Radon  Transform  (GRT)  migration  of  data  from  a  dense  array  of  broadband 
seismometers  running  west  to  east  across  central  Washington  State.  We  find  a 
low  velocity  dipping  layer  consistent  with  the  descending  oceanic  crust,  as  well 
as  a  roughly  parallel  transition  to  deeper  high  velocities  representing  the  oceanic 
Moho.  We  clearly  image  a  positive  velocity  gradient  at  a  depth  of  approximately 
40  km  beneath  the  volcanic  arc  consistent  with  the  continental  Moho.  A  fading 
velocity  contrast  within  the  oceanic  crust  consistent  with  dehydration  reactions 
starts  at  a  depth  of  ~40  km,  but  low  velocities  within  the  crust  persist  to  depths  of 
at  least  80  km  and  perhaps  100  km.  A  low  velocity  anomaly  in  the  shallow 
mantle  wedge  suggests  extensive  serpentinization. 
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We  are  able  to  compare  the  results  of  this  study  with  the  previous  GRT 
migration  studies  to  the  north  and  south.  While  the  images  are  similar  in  many 
ways,  some  differences  are  quite  striking.  Specifically,  the  fading  velocity 
contrast  normally  associated  with  dehydration  reactions  along  the  path  towards 
eclogitization  is  less  pronounced  in  the  central  Washington  profile,  despite 
presumably  similar  thermal  profiles.  Instead,  much  of  the  low  velocity  signature 
is  retained  to  depths  of  more  than  80  km.  Additionally,  we  find  a  low  velocity 
region  above  the  slab  at  this  depth,  consistent  with  a  significant  fluid  or  melt 
pocket  that  is  also  not  evident  in  the  other  profiles.  We  argue  that  for  the  central 
Washington  line,  hydrated  minerals  in  the  subducting  upper  mantle  release 
fluids  that  accelerate  the  eclogitization  of  metastable  gabbro  in  the  lower  crust, 
an  interpretation  that  also  might  help  explain  the  presence  of  volcanoes  such  as 
Mt.  Rainier  and  Mt.  Helens  in  the  Cascades  forearc. 


1.  Introduction 

Slab  derived  fluids  released  in  shallow  to  moderate  (40-130  km)  subduction 
zone  settings  play  a  critical  role  in  the  subduction  process.  The  transformation  of 
hydrated  metabasalts  to  eclogite  in  the  subducting  oceanic  crust  has  been  linked 
to  intraslab  earthquakes  (Hacker  et  al.,  2003b),  and  the  accompanying  12-15% 
increase  in  density  may  be  associated  with  changes  in  slab  dip  (Klemd  et  al., 

2011,  Bostock  et  al.,  2002)  and  increasing  stresses  within  the  slab  (Kirby  et  al., 
1996).  A  significant  portion  of  the  fluids  released  during  these  dehydration 
reactions  typically  make  their  way  into  the  mantle  wedge,  lowering  the  melting 
temperature  of  the  mantle  rock  and  potentially  resulting  in  partial  melting  and 
arc  volcanism.  The  released  fluids  also  drive  hydration  reactions  such  as 
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serpentinization  in  parts  of  the  mantle  wedge  that  are  sufficiently  cool.  This 
reduces  the  mechanical  strength  of  the  boundary  between  the  slab  and  wedge 
(Hyndman  and  Peacock,  2003),  and  can  result  in  the  creation  of  an  aseismic 
region  along  the  weakened  interface  as  the  seaward  corner  of  the  mantle  wedge 
becomes  decoupled  from  the  descending  slab  (Wada  et  al.,  2008,  Furukawa, 

1993). 

It  is  generally  recognized  that  earthquakes  occurring  at  depths  greater  than  35 
km  in  subduction  settings  are  probably  associated  with  a  process  in  which 
internal  pore  pressure  generated  by  fluids  released  during  dehydration  reactions 
reduces  the  effective  normal  stresses  along  pre-existing  faults  (Kirby  and  Wang, 
2002,  Hacker  et  al.,  2003b).  In  Cascadia,  fluid  release  from  hydrated  metabasalts 
during  eclogitization  reactions  occurs  at  depths  near  -35-45  km  (Hacker  et  al., 
2003a,  van  Keken  et  al.,  2011,  Bostock,  2012). 

Dehydration  of  hydrated  subducted  oceanic  mantle  potentially  provides  an 
additional  source  of  fluid  release  at  moderate  depths,  generally  occurring  near 
the  600°  C  isotherm  for  antigorite  and  the  800°  isotherm  for  chlorite  (Schmidt  and 
Poli,  1998,  Ulmer  and  Trommsdorff,  1995)  corresponding  to  depths  of  60-70  km 
for  the  former  and  85-95  km  for  the  latter  in  Cascadia  (van  Keken  et  al.,  2011). 
Serpentine  rocks  can  form  at  the  spreading  ridge  for  slow-spreading  and 
intermediate-spreading  centers,  and  has  been  associated  with  the  presence  of  a 
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well-developed  axial  valley  (Cannat,  1993,  Carbotte  et  al.,  2006).  Additionally, 
both  the  crust  and  upper  mantle  of  the  incoming  oceanic  plate  can  become 
hydrated  during  faulting  associated  with  bending  near  the  trench  (Ranero  et  al., 
2005,  Key  et  al.,  2012),  or  during  the  reactivation  of  spreading  center  derived 
faults  (Nedimovic  et  al.,  2008). 

In  this  study,  we  employ  a  2-D  GRT  migration  of  teleseismic  data  collected 
from  a  dense  array  of  41  seismometers  over  a  period  of  ~27  months.  We  show 
high  resolution  images  of  the  subsurface  below  central  Washington,  providing 
improved  constraints  on  the  dehydration  and  hydration  reactions  taking  place 
within  the  Cascadia  subduction  zone. 

2.  Geological  and  Geophysical  Background 

2.1  The  Cascadia  Subduction  System 

The  tectonics  of  the  Pacific  Northwest  are  dominated  by  the  Cascadia 
subduction  system.  Beneath  central  Washington,  the  Juan  de  Fuca  plate  is 
subducting  at  a  rate  of  approximately  35  mm/yr  (Miller  et  al.,  2001).  The  dip  of 
the  subducting  plate  is  between  5-15°  near  the  coast,  steepening  to  ~  40°  after 
reaching  a  depth  of  70-80  km  (McCrory  et  al.,  2006,  Roth  et  al.,  2008). 
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While  smaller  in  magnitude  than  megathrust  events,  the  impact  of  intraslab 
earthquakes  should  not  be  underestimated.  Since  1945,  there  have  been  several 
events  with  magnitudes  >5.5,  including  the  6.8  Nisqually  earthquake  in  2001  that 
resulted  in  an  estimated  monetary  loss  of  more  than  USD  two  billion  (Kirby  and 
Wang,  2002).  Most  of  these  events  have  epicenters  located  between  35  and  70  km 
depth,  with  some  small  events  reaching  depths  up  to  100  km  (Preston  et  al., 
2003).  Additionally,  the  Puget  lowland  region  of  central  and  northern 
Washington,  which  is  a  focus  of  this  study,  experiences  a  much  higher  incidence 
of  these  earthquakes  than  neighboring  regions  to  the  north  and  south. 

The  subduction  zone  also  supports  a  chain  of  active  volcanoes,  extending 
from  Mt.  Meager  in  British  Columbia  to  the  Lassen  region  in  California.  While 
the  majority  of  the  Cascade  volcanoes  lie  along  the  Quaternary  axis,  including 
those  relevant  to  the  previous  studies  in  British  Columbia  and  Oregon,  two  large 
volcanoes  and  two  smaller  vent  fields  lie  well  into  the  forearc,  as  much  as  50  km 
to  the  west  of  the  axis.  One  of  these  volcanoes,  Mt.  St.  Helens,  erupted 
catastrophically  on  May  18,  1980,  resulting  in  57  deaths  and  damages  of  more 
than  USD  one  billion.  The  other,  Mt.  Rainier,  is  a  primary  target  of  this  study. 

2.2  Previous  geophysical  studies 
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The  potential  hazards  and  proximity  to  major  urban  centers  has  predictably 
led  to  extensive  geophysical  efforts  directed  towards  a  better  understanding  of 
the  regional  substructure  and  its  dynamics.  Regional  P-wave  tomographic 
studies  have  mapped  the  extent  of  the  Juan  de  Fuca  plate  (Roth  et  alv  2008),  and 
the  velocity  structure  of  the  plate,  megathrust,  forearc  crust,  and  upper  mantle 
(Ramachandran  et  al.,  2006),  arguing  that  the  intraslab  earthquakes  occurring  at 
depths  of  40-55  km  depth  are  near  the  oceanic  Moho.  The  descending  slab  has 
been  imaged  to  depths  of  350  km  below  Washington  (Roth  et  al.,  2008,  Schmandt 
and  Humphreys,  2010)  but  shows  a  weak  or  even  non-existent  signature  below 
160  km  under  much  of  Oregon  (Schmandt  and  Humphreys,  2010). 

Newer  techniques  such  as  ambient  noise  and  multiple  plane-wave  surface 
wave  tomography  have  also  been  used  in  isolation  (Yang  et  al.,  2008)  or  in 
conjunction  with  receiver  functions  (Gao  et  al.,  2011)  to  map  the  S-velocity 
structure  of  the  upper  75-150  km,  with  particular  attention  to  mapping  the 
continental  Moho  in  the  region.  Receiver  functions  have  also  been  utilized 
(Eagar  et  al.,  2010)  to  map  the  deeper  410  and  660  km  discontinuities  throughout 
the  Pacific  Northwest. 

Single  station  methods  have  been  combined  with  GPS  measurements  to 
investigate  episodic  tremor  and  slip  (ETS)  (Rogers  and  Dragert,  2003,  Brudzinski 
and  Allen,  2007,  2010,  Audet  et  al.,  2010),  establishing  that  the  region  can  be 
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subdivided  into  three  or  more  zones  based  on  different  recurrence  intervals  of 


these  phenomena,  with  implications  for  stress  buildup  and  megathrust 
earthquake  periodicity. 

Array  migration  methods  have  been  previously  employed  in  the  region  as 
well,  including  generalized  radon  transform  (GRT)  migration  surveys  to  the 
north  of  the  current  study  across  Vancouver  Island  and  into  British  Columbia 
(Nicholson  et  al.,  2005),  and  to  the  south  through  central  Oregon  (Rondenay  et 
al.,  2001).  A  preliminary  GRT  study  using  a  subset  of  the  data  used  for  the 
current  work  was  also  published  (Abers  et  al.,  2009),  examining  the  relationship 
between  ETS  and  a  sharp  velocity  transition  interpreted  to  be  a  layer  of  weak 
sediment  or  overpressured  fault  zone. 

While  tomographic  methods  are  sensitive  to  volumetric  changes  in  material 
properties,  the  strength  of  the  GRT  migration  method  is  its  sensitivity  to  abrupt 
velocity  transitions  (Rondenay,  2009),  allowing  it  to  be  used  effectively  to 
constrain  the  boundaries  along  which  such  transitions  occur,  such  as  at  the  top 
of  a  subducting  slab,  or  at  the  oceanic  or  continental  Moho.  Because  dehydration 
and  hydration  reactions  such  as  eclogitization,  serpentinization,  and  serpentinite 
dehydration  also  have  characteristic  velocity  changes  associated  with  them 
(Hacker  et  al.,  2003a),  the  boundaries  of  the  regions  within  which  these  reactions 
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occur  can  also  be  constrained  (Bostock  et  al.,  2002),  providing  insight  into  the 
fluid  processes  taking  place  within  the  subduction  setting. 

3.  Method 

3.1  Generalized  Radon  Transform 

In  this  study,  we  use  a  2-D  generalized  radon  transform  (GRT)  migration 
approach  (Miller  et  al.,  1987,  Bostock  et  al.,  2001)  to  produce  a  series  of  high- 
resolution  seismic  images  across  the  Cascadia  subduction  zone.  The  specifics  of 
this  method  have  been  described  in  detail  in  other  papers  (see,  e.g.  Bostock  et  al., 
2001,  Shragge  et  al.,  2001,  Rondenay  et  al.,  2001,  Rondenay,  2009),  as  have  the 
associated  assumptions,  resolution  and  image  robustness  (Rondenay  et  al.,  2005, 
Rondenay  2009),  so  we  present  here  only  a  brief  review  of  the  main 
characteristics. 

The  principle  behind  migration  is  rooted  in  the  idea  that  any  point  within  the 
Earth  at  which  scattering  occurs  can  be  identified  by  back-propagation  to  depth 
of  the  scattered  wavefield  recorded  by  a  dense  array  of  receivers  at  the  surface, 
such  that  the  transposed  wavefield  will  highlight  the  sources  of  the  seismic 
scattering  in  the  subsurface.  The  Green's  function  of  the  scattered  wavefield  is 
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expressed  as  a  volume  integral,  then  simplified  by  assuming  single  scattering 
(Born  approximation)  and  the  high  frequency  ray  approximation  (Bostock  and 
Rondenay,  1999).  The  method  combines  the  stacking  of  the  scattered  wavefields 
along  diffraction  hyperbolae  with  an  inversion/backprojection  operator  based  on 
analogy  with  the  generalized  Radon  Transform  (Miller  et  al.,  1987,  Bostock  et  al., 
2001). 

In  practice,  the  2-D  GRT  inversion  allows  for  the  simultaneous  treatment  of 
five  independent  modes  of  scattering,  the  forward  scattered  PS  mode,  and  four 
backscattered  modes  -  PP,  PS,  SSh,  and  SSv.  The  travel  paths  for  each  of  these 
modes  are  depicted  in  Fig.  1. 

3.2  Preprocessing 

Pre-processing  of  the  data  was  accomplished  through  the  application  of  the 
following  series  of  steps  (see,  e.g.,  Rondenay  et  al.,  2005). 

1)  Rotation  and  transformation  of  the  waveforms  into  radial,  transverse,  and 
vertical  components 

2)  Application  of  the  inverse  free  surface  transfer  matrix  to  rotate  the  wavefield 
into  the  polarization  directions  of  the  incident  P-wave  and  the  corresponding  S 
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waves  (Sv  and  Sh),  while  suppressing  the  downgoing  free-surface  reflections 


(Kennett,  1991). 

3)  Alignment  of  the  P  traces  using  a  multi-channel  cross  correlation  (VanDecar 
and  Crosson,  1990)  and  manual  filtering  to  remove  any  excessive  low  frequency 
noise. 

4)  Application  of  eigenimage  decomposition  (Bostock  and  Sacchi,  1997,  Freire 
and  Ulrych,  1998)  to  identify  the  correlated  portion  of  the  signal,  making  the 
assumption  that  the  correlated  signal  is  a  good  estimate  of  the  incident 
wavefield.  The  uncorrelated  portion  of  the  signal  is  taken  to  be  the  scattered 
wavefield  in  the  direction  of  incident  polarization. 

5)  Deconvolution  of  the  estimated  incident  wavefield  from  the  components  of 
the  scattered  wavefield  (Berkhout,  1977,  Rondenay,  2005)  to  obtain  the 
normalized  scattered  wavefield  from  an  impulse  response.  A  station  specific 
damping  factor  (Pearce  et  al.,  2012)  was  used  so  that  damping  of  relatively  noisy 
stations  would  not  result  in  a  loss  of  information  at  other  stations 

6)  Recasting  of  the  scattered  wave  impulse  responses  into  the  reference  frame 
dictated  by  the  2-D  geometry  of  the  target  structure  and  convolution  of  the 
resulting  data  sections  with  a  phase  shifting  filter  to  recover  material  property 
perturbations  (Bostick  et  al.,  2001,  Rondenay  et  al.,  2005). 
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3.3  Modes 


The  migration  of  the  scattered  wavefield  is  performed  for  a  number  of 
different  modes  as  described  in  the  GRT  section  (3.1)  and  Fig.  1.  The  image  for 
an  individual  mode  is  obtained  by  inverting  scattered-wave  travel  times 
associated  with  that  mode,  subject  to  assumptions  made  about  the  background 
model  (1-D  isotropic  background  velocity  model,  2-D  geometry  of  scatterers). 
While  single-mode  images  can  provide  important  information  about  subsurface 
features,  it  is  important  to  be  cognizant  of  possible  cross-contamination  by  other 
modes  of  scattering.  This  occurs  because  the  peaks  in  the  coda  associated  with 
other  modes  are  interpreted  as  earlier  or  later  signals,  and  manifest  themselves  as 
shallower  or  deeper  versions  of  the  actual  features. 

The  composite  images  are  generated  by  combining  multiple  modes  into  a 
single  image.  Here,  the  problem  of  cross-contamination  is  reduced,  because  the 
signal  corresponding  to  real  structure  will  sum  up  constructively,  whereas  the 
multiples  and  signals  from  other  modes  will  not,  leading  to  a  real  structure 
signal  that  is  more  focused  and  exhibits  higher  amplitudes  (Fig.  1).  Conversely, 
as  the  modes  do  not  all  exhibit  identical  resolution,  composite  images  sometimes 
suffer  from  the  inclusion  of  modes  that  are  more  poorly  resolved. 
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4.  Data 


Our  high-resolution  seismic  images  were  generated  using  data  from  a  dense 
array  of  41  three-component  broadband  seismometer  stations,  deployed  along  a 
roughly  west  to  east  profile  across  Washington  state  (Fig.  2c).  This  array  was 
operational  from  July  2006  to  September  2008,  although  not  every  station  was 
recording  for  the  entire  time.  Thirty -five  of  the  stations  were  part  of  the  Cascadia 
Array  for  Earthscope  (CAFE)  experiment,  five  were  part  of  Earthscope's 
Transportable  Array  (TA),  and  the  remaining  station  was  operated  by  the  Pacific 
Northwest  Seismic  Network  (PNSN).  The  raw  data  are  now  publicly  available 
from  the  archive  for  seismic  data  maintained  by  the  Incorporated  Research 
Institutions  for  Seismology  (IRIS). 

The  stations  were  deployed  along  a  dense,  narrow  band  in  the  west  to  prevent 
image  aliasing  in  the  part  of  the  system  where  the  slab  is  shallowest  (Rondenay 
et  al.,  2005,  Suckale  et  al.,  2009).  As  that  constraint  could  be  relaxed  where  the 
slab  is  deeper,  the  array  was  allowed  to  fan  out  in  the  east  to  better  accommodate 
for  other  seismic  methods  such  as  tomography.  For  the  purpose  of  2-D  GRT 
imaging,  the  stations  were  projected  along  a  line  oriented  100°  from  north  -  i.e., 
generally  perpendicular  to  the  geological  strike  of  the  subducting  slab  at  this 
location  (Abers  et  al.,  2009)  (Fig.  2a). 
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We  selected  events  with  a  minimum  magnitude  of  5.7  that  exhibited  signal  to 
noise  ratio  sufficiently  high  for  successful  alignment  of  the  incident  wavefield. 
Also  required  was  an  epicentral  distance  greater  than  30°  to  prevent  triplication 
from  the  mantle  transition  zone  (Rondenay  et  al.,  2001)  and  permit  an  incident 
plane  wave  assumption  (Bostock  et  al.,  2001)  and  less  than  98°  to  rule  out  core 
diffracted  phases.  Records  from  the  events  that  were  stable  during 
deconvolution  were  migrated,  and  only  those  reasonably  free  of  ringing  and 
exhibiting  some  sense  of  structure  were  used  for  the  composite  images.  We 
generated  migrated  images  for  individual  events  using  different  combination  of 
modes  as  well,  in  order  to  assess  what  each  mode  was  contributing  to  the  image 
for  a  given  event.  We  removed  modes  that  had  very  poor  signal  to  noise  or 
demonstrated  excessive  ringing.  In  all,  data  from  sixty -three  events  covering  a 
fairly  comprehensive  range  of  back-azimuth  were  used  to  generate  our  primary 
composite  images  (Fig.  2b).  A  complete  list  of  the  events  used  for  these 
inversions  along  with  a  breakout  of  the  modes  used  for  each  event  can  be  found 
in  the  supplemental  material  (Supp.  T.  1). 

5.  Results 
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The  2-D  GRT  images  are  presented  in  figure  3.  The  profiles  show  velocity 
perturbations  relative  to  a  background  model  as  a  function  of  horizontal  distance 
(hereafter  the  x  coordinate  axis)  and  depth,  with  blue  representing  an  increase  in 
velocity  and  red  a  decrease.  The  GRT  method  has  sensitivity  to  abrupt  changes 
in  velocity,  which  can  be  identified  in  the  images  as  transitions  from  red  to  blue 
(slow  to  fast)  or  blue  to  red  (fast  to  slow).  A  slow-to-fast  transition  with 
increasing  depth  is  referred  to  as  a  positive  discontinuity,  and  the  opposite  (fast- 
to-slow),  a  negative  discontinuity. 

The  single  mode  image  in  figure  3a  (P- velocity  or  ba/a)  was  produced  using 
the  backscattered  PP  mode.  The  composite  image  in  Fig.  3b  (S-velocity  or  5|3/(3) 
was  produced  by  combining  the  inversion  results  for  the  backscattered  Ps,  SSv, 
and  SSh  modes  along  with  the  results  for  the  forward  scattered  Ps  mode  (Fig  3c-f 
respectively,  c.f..  Fig  1  for  definition  of  modes).  The  inclusion  of  the 
backscattered  modes  greatly  improves  resolution  (Rondenay  et  al.,  2009),  and  we 
can  use  the  fact  that  each  modal  image  focuses  data  that  is  independent  from  that 
focused  by  other  modes  to  establish  an  effective  robustness  test.  Accordingly,  we 
consider  a  feature  to  be  robust  only  when  appears  not  only  in  the  composite 
image,  but  in  multiple  individual  modal  images. 

The  most  prominent  feature  in  the  composite  5|3/|3  image  (Fig.  3b)  is  a  low 
velocity  (red)  layer  (LVL)  that  dips  towards  the  east.  On  the  western  (left)  edge 
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of  the  image,  the  top  of  this  layer  is  marked  by  a  negative  velocity  discontinuity 
at  a  depth  of  22  km,  and  the  bottom  by  a  positive  discontinuity  at  a  depth  of  30 
km.  Following  previous  GRT  migration  work  conducted  in  subduction  settings, 
we  interpret  this  LVL  to  represent  the  subducting  oceanic  crust  (Rondenay  et  al., 
2001,  2005,  Nicholson  et  al.,  2005,  Abers  et  al.,  2009). 

The  positive  discontinuity  indicates  the  location  of  the  oceanic  Moho,  from 
which  we  can  calculate  a  dip  angle  for  the  descending  slab  of  ~8°  between  ~32 
and  ~42  km  depth  (x  =  0  to  ~65  km),  which  steepens  to  -20°  between  ~42  and  ~90 
km  depth  (x  =  ~65  to  -180  km).  Beyond  this  depth  the  discontinuity  can  no 
longer  be  traced.  These  values  are  consistent  with  previous  estimates  of  slab  dip 
based  of  Wadati-Benioff  seismicity  in  the  region  (McCrory  et  al.,  2006). 

The  LVL  maintains  an  apparent  thickness  of  -7-8  km  to  a  depth  of  -35  km 
(x=60  km)  where  it  appears  to  thicken  considerably.  Similar  occurrences  have 
been  identified  in  previous  work  (Rondenay  et  al.,  2001,  Bostock  et  al.,  2002, 
Nicholson  et  al.,  2005),  where  the  thickening  was  interpreted  as  serpentinization 
of  the  cold  nose  of  the  mantle  wedge.  In  the  previous  studies,  the  apparent 
thickening  was  accompanied  by  the  near  disappearance  of  the  low  velocity  crust, 
interpreted  as  the  onset  of  eclogitization  reactions  within  the  crust.  While  a 
velocity  increase  occurs  here  for  central  Washington  as  well,  the  change  is  much 
more  subdued,  reducing  the  magnitude  of  the  velocity  perturbation  by  less  than 
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a  factor  of  two  (i.e.  from  d[3/(3=~-.15  at  x<90  km  to  d(3/(3= — .10  at  x>90  km).  The 
diminished  LVL  continues  at  a  steeper  dip  until  it  becomes  indistinguishable 
from  the  background  velocities  at  a  horizontal  distance  of  x=150-200  km  and  a 
depth  of  80-95  km.  Directly  above  where  the  LVL  fades  away  completely,  we 
find  another  low  velocity  feature,  extending  between  depths  of  65-85  km  and 
from  x=150-220.  Like  the  persistence  of  the  LVL  to  depth,  this  feature  appears  to 
be  isolated  to  this  part  of  the  Cascadia  subduction  system 

An  additional  robust  feature  observed  in  Fig.  3b  is  a  positive  velocity 
discontinuity  at  an  average  depth  of  ~42  km,  which  extends  from  x=180  to  x=235. 
Again  following  previous  work  (Rondenay  et  al.,  2001,  Nicholson  et  al.,  2005),  we 
interpret  this  feature  as  the  Moho  of  the  overriding  continental  plate.  The 
velocity  contrast  is  not  as  strong  here  as  it  was  in  those  studies  and  in  fact  was 
not  evident  at  all  in  previous  work  using  only  part  of  the  current  data  set  (Abers 
et  al.,  2009).  The  discontinuity  appears  to  exhibit  some  topography,  including  a 
concave  up  segment  between  x=l 80-200  followed  by  a  concave  down  segment 
between  x=200-220. 

The  5  (3/(3  images  are  generally  considered  to  be  more  robust  than  the  ba/a 
images  due  both  to  a  more  complete  separation  of  the  S-scattered  waves  from  the 
incident  wavefield  and  to  the  availability  of  multiple  modes  which  improve 
volume  and  dip  resolution  (Rondenay,  2009).  An  examination  of  the  ba/a 
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images  is  still  warranted  as  it  can  provide  additional  and  independent 
information.  In  particular,  they  can  be  used  to  establish  the  robustness  of  certain 
structures,  since  multiple  contamination  is  less  prevalent  in  the  5a/a  images 
(Pearce  et  al.,  2012).  Moreover,  differences  between  the  5|3/(3  and  ba/a  images 
can  highlight  regions  in  the  subsurface  where  the  background  velocity  model 
breaks  down.  If  the  background  model  is  inaccurate  with  respect  to  the  Vp/Vs 
ratio,  the  features  will  be  mapped  at  the  different  depths  on  each  image. 

Our  primary  ba/a  image  (Fig.  3a)  exhibits  many  of  the  same  features  as  the 
5|3/(3  image.  The  positive  velocity  transition  at  ~42  km  in  the  east  is  present, 
extending  from  x=170  km  to  210  km  along  the  horizontal  axis.  The  topography 
of  this  interface  noted  in  the  5|3/|3  image  is  arguably  present  here  as  well,  though 
not  as  well  defined.  The  dipping  low  velocity  layer  in  the  west  is  present,  similar 
in  location,  thickness,  and  dip  to  that  observed  in  the  5 [3/p  image.  The  extension 
of  the  low  velocity  layer  becomes  difficult  to  trace  beyond  a  depth  of  ~65  km, 
likely  due  to  the  limited  resolution  of  the  ba/a  image,  which  becomes 
increasingly  limited  with  depth  (Rondenay  et  al.,  2005). 

6.  Discussion 

6.1  Central  Washington 
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6.1.1  The  continental  Moho 


We  interpret  the  positive  velocity  discontinuity  at  a  depth  of  ~42  km  in  the 
eastern  part  of  the  image  as  the  continental  Moho,  consistent  with  previous  GRT 
migration  studies  (Rondenay  et  al.,  2001,  Bostock  et  al.,  2002,  Nicholson  et  al., 
2005,  Rondenay  et  al.,  2010).  In  a  preliminary  study  using  a  subset  of  this  data 
(Abers  et  al.,  2009),  the  continental  Moho  was  not  well  imaged,  probably  due  to 
diminished  station  coverage,  as  data  from  two  critical  stations  were  not 
recovered  in  time  for  that  study  due  to  weather  conditions  in  the  mountains. 

The  topography  that  appears  along  the  continental  Moho  is  consistent  with 
isostatic  compensation  in  response  to  the  influence  of  the  nearby  Cascade  Range, 
which  is  located  directly  above  the  concave  up  feature  in  the  discontinuity.  The 
6  km  deflection  in  the  Moho  would  be  consistent  with  a  1.1  km  rise  in  the  surface 
topography  given  a  mantle  density  of  3300  kg/m3  and  a  crustal  density  of  2800 
kg/m3.  This  1.1  km  rise  is  consistent  with  the  local  elevation,  and  is  roughly  a 
quarter  of  Mt.  Rainier' s  peak  height  of  4.17  km. 

6.1.2  The  subduction  system 
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We  interpret  the  dipping  low  velocity  layer  extending  from  the  western  edge 
of  the  image  as  the  subducting  oceanic  crust.  The  layer  has  a  thickness  of  ~8  km 
and  dips  at  slightly  less  than  10°,  consistent  with  what  we  would  expect  for  the 
oceanic  crust  at  this  location.  The  strong  downward  negative  velocity  gradient  at 
the  top  of  the  layer  has  been  interpreted  to  be  either  a  weak  layer  of  sediment  or 
fluids  in  pressurized  channels  (Abers  et  al.,  2009). 

The  positive  velocity  contrast  at  the  bottom  of  the  layer  has  typically  been 
interpreted  to  represent  the  oceanic  Moho  (Abers  et  al.,  2009),  which  can  be  seen 
reasonably  clearly  to  a  depth  of  more  than  90  km.  This  velocity  contrast  is 
reduced  in  multiple  stages,  first  at  a  depth  of  ~40  km,  and  again  at  a  depth  of  ~80 
km,  corresponding  to  increases  in  the  velocity  of  the  subducting  crust.  The 
velocity  increase  at  a  depth  of  ~40  km  is  interpreted  to  reflect  a  material 
transformation  in  the  upper  crust  from  hydrated  metabasalts  to  zoisite- 
amphibole  eclogite  (Helffrich,  1996,  Bostock  et  al.,  2002,  Hacker  et  al.,  2003a), 
along  with  a  reduction  in  fluid  pore  pressure  corresponding  to  a  rupturing  of  the 
plate  boundary  seal  (Audet  et  al.,  2010,  Bostock,  2012).  For  warm  subduction 
zones,  the  eclogitization  of  hydrated  metabasalts  is  largely  controlled  by  pressure 
(Schmidt  and  Poli,  1998,  Hacker  et  al.,  2003a,  Rondenay  et  al.,  2008)  and  is 
expected  to  occur  at  depths  of  -40-55  km  in  Cascadia  (van  Keken  et  al.,  2011). 
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The  fluids  released  from  this  eclogitization  reaction  migrate  upward  into  the 
mantle  wedge  (Peacock,  1990),  where  they  are  taken  up  by  hydration  reactions 
resulting  in  the  generation  of  several  hydrous  minerals  such  as  chlorite, 
antigorite  serpentine,  talc,  and  brucite  (Hyndman  and  Peacock,  2003,  Ulmer  and 
Trommsdorf,  1995).  Like  eclogitization,  serpentinization  of  mantle  material  has 
seismological  consequences,  including  a  significantly  lower  velocity  and  an 
increase  in  Poisson's  ratio  (Bostock  et  al.,  2002,  Hyndman  and  Peacock,  2003). 
We  interpret  the  low  velocity  region  overlying  the  40  km  deep  crust  to  be  the 
serpentinized  mantle  wedge. 

Our  interpretation  of  the  downgoing  oceanic  crust  and  Moho,  eclogitization 
reaction  within  the  oceanic  crust,  and  serpentinized  corner  of  the  mantle  wedge 
are  consistent  with  identifications  made  during  previous  GRT  migration  studies 
in  the  Cascadia  region  (Rondenay  et  al.,  2001,  Nicholson  et  al.,  2005,  Abers  et  al., 
2009).  However,  there  are  some  important  differences  between  the  images  from 
the  work  along  other  transects  and  the  one  from  the  current  study  that  must  be 
addressed. 

The  first  of  these  differences  relates  to  the  retention  of  a  significant  low 
velocity  signature  in  the  subducting  oceanic  crust  well  beyond  the  eclogitization 
of  the  hydrated  upper  crust.  This  transition  to  zoisite  amphibole  eclogite  is 
expected  to  raise  the  P-wave  velocity  of  the  crust  from  -15%  slower  than  the 
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underlying  mantle  to  ~4%  slower  at  depths  starting  at  ~  35  km,  with  a  further 
transition  at  -7 0  km  to  anhydrous  eclogite  which  is  expected  to  render  the  crustal 
velocity  signature  indistinguishable  from  the  surrounding  material  (Hacker  et 
al.,  2003a).  Our  model  shows  a  more  subdued  velocity  change  at  ~35  km, 
retaining  a  clear  contrast  at  the  bottom  of  the  layer  to  depths  of  80-95  km.  This 
contrast  can  be  explained  by  the  presence  of  a  lower  crustal  layer  of  metastable 
gabbro  (Hacker  et  al.,  2003a,  Rondenay  et  al.,  2008),  which  could  be  as  much  as  5 
km  thick  (Hacker,  2008,  van  Keken  et  al.,  2011)  and  with  a  velocity  as  much  as 
12%  slower  than  the  underlying  harzburgite  (Hacker  et  al.,  2003a). 

The  coarse-grained  metastable  gabbro  lower  crustal  layer  can  transform 
directly  to  eclogite  at  pressures  and  temperatures  that  correspond  to  depths  of 
-80-90  km  for  Cascadia  (John  and  Schenk,  2003,  Hacker  et  al.,  2003b),  but  this 
reaction  tends  to  be  extremely  sluggish  in  the  absence  of  free  water  (Hacker  et  al., 
2008).  A  key  catalyst  is  the  introduction  of  water  to  pathways  within  the  lower 
crust,  therefore  it  is  expected  that  dehydration  of  serpentine  or  chlorite  in  the 
harzburgite  of  the  underlying  mantle  would  accelerate  this  process  (Schmidt  and 
Poli,  1998,  John  and  Schenk,  2003).  The  released  fluids  would  access  the  gabbroic 
mass  through  channels  created  by  volume  reduction  associated  with  the 
initiation  of  prograde  metamorphism  or  alternately  through  stress  related 
shearing  (Bruhn  et  al.,  2000,  John  and  Schenk,  2003).  Once  initiated,  the 
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transition  from  gabbro  to  eclogite  furthers  the  volume  reduction  of  the  host  rock, 
increasing  the  permeability  and  thereby  enabling  the  infiltration  of  any 
remaining  fluid  (Aharonov  et  al.,  1997,  John  and  Schenk,  2003). 

The  second  notable  feature  in  the  image  that  must  be  explained  is  the  low 
velocity  region  that  extends  from  x=170-210  and  from  depths  of  60-80  km,  which 
is  consistent  with  a  fluid  or  melt  phase  (Rondenay  et  al.,  2010).  The  base  of  the 
feature  appears  to  be  roughly  coincident  with  the  top  of  the  slab  from  depths  of 
75-90  km  (mantle  depths  of  80-95  km),  the  same  range  over  which  the  low 
velocity  signature  of  the  subducting  crust  disappears.  While  serpentinized 
peridotite  in  the  upper  mantle  is  only  stable  to  depths  of  -60-70  km  in  Cascadia 
(Hacker  et  al.,  2003a,  van  Keken  et  al.,  2011),  chlorite  bearing  harzburgite  is  stable 
to  temperatures  approaching  800  C  (Schmidt  and  Poli,  1998,  Hacker  et  al.,  2003a), 
corresponding  to  depths  of  -80-95  km  in  Cascadia  (Hacker  et  al.,  2003a)  with 
complete  dehydration  of  the  uppermost  mantle  by  115  km  depth  (van  Keken  et 
al.,  2011).  This  suggests  that  chlorite  rather  than  antigorite  serpentine  would 
have  to  be  present  in  the  upper  mantle. 

The  thermal  profile  in  Cascadia  does  not  allow  for  extensive  hydration  of  the 
oceanic  upper  mantle,  but  evidence  of  hydration  has  been  found  within  the 
uppermost  few  km  (Kao  et  al.,  2008,  Nedimovic  et  al.,  2009).  The  velocity  of 
chlorite  harzburgite  is  roughly  midway  between  that  of  metastable  gabbro  and 
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spinel  harzburgite  (Hacker  et  al.,  2003a),  so  a  relatively  thin  and  localized 
hydrated  upper  mantle  would  blend  into  the  low  velocity  crustal  signature.  This 
would  imply  that  the  extended  low  velocity  layer  imaged  to  depths  of  80-95  km 
is  a  combination  of  a  ~4  km  thick  metastable  gabbro  layer  underlain  by  layer  of 
chlorite  harzburgite,  also  a  few  km  thick. 

Figure  5a  shows  the  location  of  earthquake  hypocenters  occurring  within  25 
km  of  the  projection  line  for  our  survey.  Most  of  the  hypocenters  that  appear  to 
be  associated  with  the  descending  slab  are  in  locations  consistent  with  what  we 
would  expect  for  the  eclogitization  of  hydrated  basalts  in  the  upper  crust  (Hacker 
et  al.,  2003b).  A  few  of  the  larger,  deeper  ones  between  x=110  and  150  km  more 
likely  represent  dehydration  of  serpentine  in  the  upper  mantle  (Hacker  et  al., 
2003a,  Kao  et  al.,  2008).  The  systematic  biases  which  affect  hypocenter  location 
tend  to  make  them  appear  as  much  as  25  km  farther  from  the  ocean  than  they 
actually  are  (Syracuse  and  Abers,  2009),  so  all  of  these  hypocenters  may  well  be 
located  somewhere  within  the  mantle  of  the  descending  slab.  There  is  a  cluster 
of  small  earthquakes  that  do  not  fit  this  profile,  at  depths  of  88-100  km  and 
x=200-225  km,  that  are  located  directly  below  the  low  velocity  fluid/melt  feature. 
As  the  transition  of  metastable  gabbro  to  eclogite  is  expected  to  be  aseismic  (John 
and  Schenk,  2003),  these  hypocenters  likely  represent  dehydration  reactions  in 
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the  oceanic  mantle,  with  depths  that  are  consistent  with  the  dehydration  of 
chlorite. 

The  introduction  of  hydration  alteration  into  the  uppermost  mantle  can  occur 
in  a  number  of  ways.  Serpentinized  mantle  rocks  have  been  found  in  oceanic 
lithosphere  at  slow  and  intermediate  spreading  ridges  (Cannat,  1993,  Carbotte  et 
al.,  2006),  but  the  narrow  axial  valley  at  the  Juan  de  Fuca  ridge  rules  out  this 
possibility  for  Cascadia  (Canales  et  al.,  2005,  Carbotte  et  al.,  2006).  Likewise, 
hydration  of  the  lower  crust  and  upper  mantle  has  been  shown  to  occur  during 
bend  related  faulting  at  the  trench  in  Central  America  (Ranero  et  al.,  2005,  Key  et 
al.,  2012),  but  the  shallow  dip  angle  of  ~7.5  degrees  at  the  trench  (McCrory  et  al., 
2006)  and  a  ~2  km  thick  sediment  load  that  provides  both  a  hydraulic  barrier  and 
thermal  blanket  (Divens  et  al.,  retrieved,  2012,  Spinelli  et  al.,  2004)  make  this 
scenario  unlikely  for  Cascadia. 

There  is  substantial  evidence  that  some  hydration  of  the  lower  crust  and 
uppermost  mantle  does  occur  in  Cascadia.  Analysis  of  the  magnitude  6.7 
Nisqually  earthquake  also  suggests  hydration  extending  as  much  as  10  km  into 
the  mantle  (Kao  et  al.,  2008),  and  seismic  reflection  studies  have  found  faulting 
and  hydration  up  to  200  km  seaward  of  the  trench,  in  zones  of  potential  plate 
weakness  associated  with  propagator  wakes  (Nedimovic  et  al.,  2009).  The 
hydration  of  these  faults  appears  to  be  largely  restricted  to  lower  crustal  levels. 
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but  in  some  cases,  the  faults  have  penetrated  the  Moho  (Nedimovic  et  al.,  2009). 
While  serpentinization  is  normally  the  most  important  hydration  mechanism 
associated  with  peridotites,  the  stability  zone  for  serpentine  is  expected  to  extend 
no  more  than  a  few  km  into  the  mantle  (Ulmer  and  Trommsdorf,  1995,  Hyndman 
and  Wang,  1995,  Nedimovic  et  al.,  2009).  The  widest  of  these  propagator  wake 
regions  in  Cascadia  is  coincident  with  the  deep  seismicity  in  the  eastern  part  of 
our  profile. 

6.2  Implications  for  regional  tectonics 

In  GRT  imaging,  the  projection  line  azimuth  must  be  chosen  to  be  as  close  to 
perpendicular  to  the  2-D  strike  of  the  structure  as  possible,  which  for  subduction 
zones  is  parallel  to  the  direction  of  slab  dip.  When  the  steepest  dip  is  in  the 
direction  of  slab  motion,  we  can  think  of  the  slab  in  cross-section  as  a  progression 
through  time,  such  that  the  condition  of  the  slab  at  some  depth  z2  is  a 
consequence  of  the  reaction  we  image  at  a  shallower  depth  zl.  When  the 
direction  of  steepest  dip  and  slab  motion  differ  significantly,  the  projection  line 
cuts  across  the  direction  of  slab  motion  and  complicates  this  perception  (Fig.  4). 

Figure  5b-g  show  our  study  area  in  map  view  broken  up  into  three  corridors, 
with  the  region  of  our  migrated  image  from  x=0  to  50  km  falling  into  the  X 
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corridor,  that  from  x=50  to  180  km  into  the  Y  corridor,  and  the  eastern  section 


from  x=180  km  into  the  Z  corridor.  The  corridors  have  azimuths  of  30  degrees 
north  of  east  (40  degrees  north  of  our  projection  line),  consistent  with  the 
direction  of  relative  plate  motion  (Fluck  et  al.,  1997,  Miller  et  al.,  2001).  The 
location  of  the  boundaries  for  these  corridors  was  directed  by  changes  in  the 
distribution  of  earthquake  epicenters,  which  are  also  shown  in  the  figures.  The 
same  biases  regarding  epicenter  location  (Syracuse  and  Abers,  2009)  apply  here 
as  they  did  for  figure  5a,  but  corrections  for  these  biases  would  not  undermine 
our  general  conclusions. 

The  segment  of  our  projection  line  that  crosses  the  Y  corridor  passes  directly 
over  the  large  concentration  of  events  associated  with  dehydration  reactions  at 
depths  of  38-70  km  (Fig.  5d-f).  These  events  densely  populate  the  Y  corridor  and 
extend  to  some  degree  northward,  but  are  almost  completely  absent  within  the  Z 
corridor.  Conversely,  all  of  the  events  that  occur  at  depths  of  greater  than  70  km 
for  this  region  fall  within  the  Z  corridor  (Fig.  5g). 

The  previous  GRT  migration  work  done  to  the  north  (Nicholson  et  al.,  2005) 
and  south  (Rondenay  et  al.,  2001)  provide  a  convenient  framework  within  which 
we  can  make  certain  assertions  regarding  regional  along-strike  variation.  Figure 
6a  depicts  a  regional  map  with  three  transects  representing  the  projection  lines 
from  three  different  GRT  migration  experiments.  The  A- A'  transect  (Fig.  6b) 
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runs  across  the  southeastern  end  of  Vancouver  Island  and  into  southwest  British 


Columbia  (Nicholson  et  al.,  2005).  The  C-C'  transect  (Fig.  6d)  is  from  a  similar 
study  conducted  across  central  Oregon  (Rondenay  et  al.,  2001).  The  B-B'  transect 
(Fig.  6c)  across  central  Washington  state  shows  the  projection  line  for  the  current 
experiment.  The  X,  Y,  and  Z  corridors  are  also  included  on  the  map,  in  addition 
to  a  Vancouver  Island  (VI)  corridor  which  includes  the  A- A'  transect,  and  an 
Oregon  (OR)  corridor  which  brackets  the  C-C'  transect.  The  T  corridor  bridges 
the  gap  between  Z  and  OR. 

The  three  migrated  images  share  a  number  of  similarities  in  terms  of  general 
structure,  but  also  show  some  striking  differences.  The  degree  of  disappearance 
of  the  low  velocity  layer  at  a  depth  of  40-45  km  is  very  pronounced  in  A- A' 

(more  than  80%  of  the  original  contrast),  a  bit  less  so  in  C-C'  (60%),  and 
significantly  more  subdued  in  B-B'  (30%),  suggesting  that  the  metastable  gabbro 
layer  identified  in  central  Washington  is  not  equally  present  along-strike.  The 
low  velocity  fluid/melt  feature  above  the  slab  from  x=160  to  x=220  is  also  clearly 
present  only  in  central  Washington,  so  the  deep,  concentrated  release  of  fluids 
from  the  upper  mantle  is  probably  a  local  feature  as  well,  a  conclusion  supported 
by  the  distribution  of  earthquake  hypocenters. 

While  several  parameters  vary  along-strike  in  Cascadia  (see  Table  2), 
including  subduction  rate  (Miller  et  al.,  2001),  plate  age  (Davis  and  Hyndman, 
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1989),  and  sedimentation  (Divins  et  al.,  retrieved  2012),  one  of  the  parameters 
that  correlate  best  with  seismic  activity  is  slab  dip  (Fig  6a).  Because  the  slab  dip 
is  shallow  for  the  VI,  X,  and  Y  corridors  to  depths  of  30  km,  the  megathrust 
seismological  locked  and  transitions  zones  extend  eastward,  the  latter  on  to  land 
(Fluck  et  al.,  1997).  The  shallow  earthquakes  (33-38  km  depth.  Fig.  5b)  near  the 
coast  correlate  very  well  with  the  megathrust  transition  zone,  and  the  vast 
majority  of  the  intraslab  activity  takes  place  along  these  corridors. 

The  deeper  earthquakes  occur  within  the  Z  corridor,  where  a  significant 
change  in  dip  is  accommodated.  From  a  depth  of  70  km  to  110km,  the  dip 
increases  within  the  corridor  as  much  as  8°  over  a  corridor-perpendicular 
distance  of  less  than  50  km  (McCrory  et  al.,  2006,  also  see  Fig.  6a  and  Table  2). 
The  relationship  between  the  earthquakes  and  the  change  in  dip  is  perhaps  less 
clear. 

One  possibility  is  that  the  change  in  dip  is  simply  a  result  of  deformation 
caused  by  the  simple  geometric  constraints  imposed  by  the  curvature  of  the 
subduction  zone  (Creager  and  Boyd,  1991).  While  stresses  associated  with  the 
changes  could  connect  otherwise  isolated  pockets  of  fluid  (Mibe  et  al.,  2003, 
Bruhn  et  al.,  2000)  and  the  valley  like  topography  formed  by  the  changes  in  dip 
could  channel  fluids  over  a  short  horizontal  profile,  the  presence  of  fluids  would 
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not  be  required.  Likewise,  this  interpretation  would  not  explain  differences  in 
the  Z  corridor  seaward  of  the  change  in  dip. 

Conversely,  the  changes  in  dip  could  be  a  result  of  the  densification  of  the  slab 
(Klemd  et  al.,  2011)  due  to  the  fluid  release  from  chlorite  harzburgite  in  the 
upper  mantle  and  subsequent  triggering  of  the  eclogitization  reaction  in  the 
lower  crust.  This  would  require  a  fundamental  difference  in  the  incoming  slab 
within  the  Z-corridor  that  could  enable  sufficient  hydration  of  the  upper  mantle, 
such  as  zones  of  weakness  associated  with  propagator  wakes  (Nedimovic  et  al., 
2009). 

The  largest  of  these  propagator  wake  zones  includes  the  segment  of  the  Z 
corridor  imaged  by  this  study  (Fig.  6a),  and  also  serves  to  delineate  three  very 
different  regions  of  the  Cascadia  subduction  zone.  To  the  north,  the  dip  of  the 
plate  is  shallow,  intraslab  seismic  activity  is  dense,  and  the  127  quaternary 
volcanic  vents  are  tightly  concentrated  around  the  five  major  volcanoes  (Mt. 
Meager  is  off  the  map  to  the  north) (Flildreth,  2007).  To  the  south,  the  plate  dips 
more  steeply,  seismic  activity  is  sparse,  and  more  than  2000  quaternary  vents 
form  a  near  continuous  array  from  Bumping  Lake  to  the  east  of  Mt.  Rainier  to 
very  near  the  California  border  (Flildreth,  2007). 

While  most  of  the  volcanism  in  the  Cascades  lies  along  the  quaternary  axis, 
there  are  four  volcanic  features  within  the  forearc;  three  of  them  (Mt.  Rainier,  Mt. 
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St.  Helens,  and  the  Portland  vent  fields)  lie  within  the  wake  propagator  zone, 
and  the  fourth  (Indian  Head  vent  field)  just  along  its  southern  edge. 

Another  possibility  is  that  the  changes  in  dip  are  connected  to  both  ideas,  with 
the  weaknesses  inherent  in  the  wake  propagator  zone  and  the  slab  densification 
selecting  the  Z  corridor  for  accommodation  of  the  geometric  constraints  imposed 
by  the  curvature  of  the  subduction  zone. 

7.  Conclusion 

In  this  study,  we  present  a  high  resolution  seismic  cross  section  of  the 
Cascadia  subduction  setting  across  central  Washington  State.  We  see  an  ~8  km 
thick  low  velocity  layer  descending  from  the  west  that  we  interpret  as  the 
oceanic  crust,  underlain  by  a  sharp  velocity  transition  that  we  take  to  be  the 
oceanic  Moho.  A  velocity  reduction  within  the  slab  crust  at  a  depth  of  -35-45  km 
is  consistent  with  the  transformation  of  hydrated  metabasalts  in  the  upper  crust 
to  eclogite,  and  the  low  velocity  nose  in  the  mantle  wedge  is  characteristic  of 
serpentinite  formed  in  the  continental  mantle  by  fluids  rising  from  the  slab 
during  the  eclogitization  process.  We  also  image  a  positive  velocity  transition  at 
a  depth  of  -42  km  in  the  east  that  we  interpret  to  represent  the  continental  Moho. 
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The  transition  of  the  low  velocity  oceanic  crust  due  to  the  eclogitization  of 
the  hydrated  basalts  is  more  subdued  than  expected,  suggesting  that  a  layer  of 
metastable  gabbro  is  present  in  the  lower  crust.  A  velocity  increase  within  the 
crust  starting  at  ~80  km  depth,  and  a  overlying  low  velocity  feature  consistent 
with  a  pocket  of  released  fluid  or  melt  suggests  the  dehydration  of  chlorite  in  the 
upper  mantle. 

The  region  imaged  by  our  profile  at  this  depth  appears  to  be  unique  for 
Cascadia,  as  it  corresponds  to  the  accommodation  of  a  significant  along-strike 
change  in  dip,  a  set  of  small  earthquakes  at  a  depth  of  88-100  km,  and  the 
landward  extension  of  a  wake  propagator  zone.  The  interplay  between  these 
features  could  have  important  implications  for  the  volcanic  and  seismic 
expressions  of  the  Cascadia  subduction  system. 
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FIGURES 


Figure  1:  A  simple  illustration  of  the  modes  used  for  the  migrated  sections  in  this 
paper.  The  incident  plane  P-wave  encounters  a  scattering  contact  in  the 
subsurface,  such  that  the  energy  of  the  transmitted  wave  is  divided  into  P  wave 
components  (solid  lines)  and  S  wave  components  (dashed  lines).  The  5a/a  (P- 
wave)  sections  show  the  backscattered  PP  mode  in  isolation,  whereas  the  5(3/ [3  (S- 
wave)  sections  show  composite  images  that  incorporate  the  forward  PS,  and 
backscattered  PS,  SSv,  and  SSh  modes.  The  SSv  and  SSh  modes  differ  only  in 
polarization,  with  the  SSv  being  the  component  of  the  wave  polarized  in  the 
vertical,  and  SSv  being  the  component  polarized  in  the  horizontal. 
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Figure  2:  Data:  The  bottom  of  this  figure  shows  the  regional  map  of  Washington 
state  for  the  CAFE  line,  with  the  location  of  the  41  stations  plotted.  The  names 
for  the  specific  station  and  their  position  along  the  selected  projection  line  of  100° 
east  of  north  is  given  in  the  above  left,  color  coded  to  denote  affiliation.  The  63 
events  used  to  generate  our  primary  migrated  sections  are  plotted  on  the 
hemispherical  map  in  the  upper  right. 
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Figure  3:  Raw  sections  for  the  CAFE  line  migration.  Figure  3a  is  the  primary 
5a/a  image  for  the  set,  and  figure  3b  is  the  primary  composite  6 (3/ (3  image, 
composed  of  the  backscattered  PS  (figure  3c),  the  backscattered  SSv  (figure  3d), 
the  backscattered  SSh  (figure  3e),  and  the  forward  scattered  PS  (figure  3f)  modes. 
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Figure  4:  A  comparison  of  the  consequences  of  differences  between  the  direction 
of  the  projection  line  (which  is  associated  with  steepest  dip)  and  the  direction  of 
plate  motion.  A  migrated  section  that  follows  a  projection  line  that  differs  from 
the  direction  of  plate  motion  is  sampling  different  source  regions  across  its 
length. 
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Figure  5:  Earthquake  distribution:  Figure  5a  shows  the  distribution  of 
earthquakes  across  the  migrated  image,  including  earthquakes  occurring 
between  1970-2012  within  25  km  to  either  side  of  the  projection  line.  The  size  of 
the  circles  corresponds  to  the  magnitude  of  the  earthquakes,  with  the  smallest 
representing  magnitudes  between  1.0  and  2.0,  and  the  largest  showing  the 
location  of  the  2001  6.8  magnitude  Nisqually  earthquake.  Figures  5b-5g  show  the 
earthquake  locations  in  map  view  for  a  selection  of  depths,  and  also  the  choice  of 
corridors  to  understand  the  distribution  of  earthquakes  in  the  context  of 
direction  of  plate  motion.  All  earthquake  hypocenter  locations  courtesy  of  IRIS 
seismiquery. 
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Figure  6:  Regional  interpretation:  Figure  6a  shows  a  map  of  the  region 
containing  three  GRT  migration  studies.  The  northernmost.  A- A',  is  from  a 
study  done  by  Nicholson  et  al.,  2005.  The  southernmost,  C-C',  is  from  a  study 
done  by  Rondenay  et  al.,  2001.  The  CAFE  line,  B-B'  reflects  the  work  done  for 
this  current  experiment.  The  map  includes  the  selection  of  corridors  for  the 
entire  region,  extending  from  the  trench,  and  a  plot  of  sediment  thickness  at  the 
trench  (Divens  et  al.,  retrieved  2010).  It  also  shows  slab  depth  contour  lines  for 
20,  50,  80,  and  110  km  depths  (McCrory  et  al.,  2006).  Earthquakes  shallower  than 
70  km  depth  are  shown  as  blue  symbols,  and  those  greater  than  70  km  are  shown 
as  green  symbols  (Iris  Seismiquery,  2012).  Red  triangles  mark  the  location  of 
volcanoes.  Wake  propagator  zones  are  included  as  shaded  regions  (Nedimovic 
et  al.,  2007).  Figure  6b-6d  show  the  migrated  sections  for  A- A',  B-B',  and  C-C' 
respectively.  HC=hyd  rated  crust,  OM=  oceanic  Moho,  SW  =  serpentinized 
wedge,  E=  eclogitization  reaction,  CM=  continental  Moho,  F=  fluid/melt  phase 
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TABLES 


Background  velocity  model 

mf 

Depth(km) 

Vp 

Vs 

Density 

0-1 

5.4 

3.09 

2.6 

1-4 

5.4 

3.09 

2.6 

4-9 

6.38 

3.65 

2.6 

9-16 

6.59 

3.77 

2.8 

16-20 

6.73 

3.85 

3.1 

20-25 

6.86 

3.92 

3.1 

25-41 

6.96 

3.98 

3.1 

41-71 

7.8 

4.46 

3.5 

71-300 

7.8 

4.46 

3.5 

Table  1:  The  background  velocity  model  used  for  the  final  migrated  section  for 
the  CAFE  line. 
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Rate  of 

subduction 

« 

Sediment 
depth  at 
trench  (b) 

Plate  age 
at  trench 

(C) 

Earthquakes 
20-30  km  depth 

w 

Earthquakes  30- 
44  km  depth 

(d) 

Earthquakes 
44-70  km 
depth 

(d) 

Earthquakes 
>70  km  depth 

(4) 

Initial 

subduction 
angle  (trench 
to  10km  depth) 

is) 

Subduction 
angle  prior  to 
eclogitization 
(10-30  km 
depth)  (e) 

Subduction  angle 
after 

eclogitization 
(50-70  km  depth) 
(e) 

VI 

36.9  mm/a 

1600-1800 

m 

(south  of 
Nitenat  fen) 

5-7  Ma 

Some 

Some 

35-38  km  depth 
(transition  zone 
near  coast) 
more  inland 

few 

1  near  border 
at  73.3  km 
depth 

7.1  B 

10.2  B 

21.8  B 

XY 

34.3  mm/a 

1800-2000 

m 

(Astoria 

Fan) 

7-8  Ma 

Some  in  X 

Many  in  Y 

Same  as  VI 

many 

none 

7.5  B 

8.5  B 

20.5  B 

Z 

34.3  mm/a 

2000  m 
(Astoria 

Fan) 

8  Ma 

Some 

Skewed  west 
vs.  XY 

Essentially 

none 

Essentially 

none 

2  at  72-76  and 

cluster  of  8  at 
>88  km  depth 

10.1  B 

8.7  B 

22.1  B 

T 

32.1  mm/a 

1500-1800 

m 

8.5  Ma 

Some 

Skewed  more 

west 

(follows  slab 
depth  45-50 
km) 

none 

Essentially 

none 

none 

11. 0B 

9.6  B 

18.0  B 

OR 

30.9  mm/a 

600-1400  m 

8.5-9  Ma 

Only  on  T 
border 

none 

Essentially 

none 

none 

10.1  B 

9.7B 

15  .8  B 

(a)  Miller  et  al.,  2001 

(b)  Divins  et  al..  retrieved  2012.  Fluck  et  al.,  1997 

(c)  Davis  and  Hyndman.  1989 

(d)  IRIS  seismiquery 

(e)  McCrory  et  aL,  2006 


Table  2:  A  selection  of  parameters  related  to  along-strike  variation  in  the 
Cascadia  subduction  system 
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Supplemental  Table  1:  This  table  shows  the  events  used  and  modes  used  for  each 
event  in  generating  the  primary  images.  The  events  are  in  the  year-Julian  date¬ 
time  format  used  for  IRIS  seismiquery.  Modes  designated  ‘Y  were  used  in  the 
composite  images,  while  modes  designated  'O'  were  not.  The  last  column  gives 
the  quadrant  for  the  backazimuth  of  a  given  event. 


EVENT 

S 

PP  PS  SSv  SSh  Quadrant 

’20062192218’ 

’0’ 

T 

1’ 

T 

1’ 

III 

’20062231430’ 

T 

’O’ 

T 

T 

T 

II 

’20062291111’ 

’0’ 

T 

T 

T 

’0’ 

IV 

’20062291520’ 

1’ 

T 

T 

T 

1’ 

IV 

’20062362150’ 

’0’ 

T 

T 

’O’ 

1’ 

IV 

’20062441018’ 

’0’ 

T 

T 

T’ 

’0’ 

III 

’20062710622’ 

’1’ 

’O’ 

’O’ 

T’ 

’O’ 

III 

’20062721308’ 

1’ 

T 

’O’ 

1’ 

1’ 

II 

’20062731750’ 

’0’ 

T 

T’ 

1’ 

T 

IV 

’20062761803’ 

’0’ 

T 

’O’ 

T 

’0’ 

III 

’20062900125’ 

’0’ 

T 

T’ 

T’ 

’0’ 

III 

’20062962117’ 

’0’ 

T 

T’ 

1’ 

1’ 

IV 

’20063191114’ 

’0’ 

T 

T 

1’ 

T 

IV 

’20063411910’ 

’1’ 

T 

T 

T 

’1’ 

IV 

’20070130423’ 

1’ 

T 

T’ 

1’ 

1’ 

IV 

’20070302137’ 

’0’ 

T 

T’ 

1’ 

T’ 

IV 

’20070770211’ 

T 

’O’ 

’O’ 

T’ 

T 

II 

’20070950356’ 

’1’ 

T 

T’ 

T 

T’ 

I 

’20071502022’ 

’0’ 

’O’ 

’O’ 

1’ 

’0’ 

IV 

’20071641929’ 

’1’ 

T 

T 

1’ 

T 

II 

’20071960927’ 

’0’ 

’O’ 

’O’ 

T 

’O’ 

III 

’20071970113’ 

’0’ 

’I’ 

T’ 

T’ 

T 

IV 

’20071971417’ 

1’ 

T 

T’ 

1’ 

1’ 

IV 

’20072272022’ 

1’ 

’O’ 

’O’ 

1’ 

’0’ 

IV 

’20072280516’ 

’1’ 

’1’ 

T 

T 

T’ 

II 

’20072280839’ 

’1’ 

1’ 

T 

T 

’I’ 

III 

’20072300252’ 

1’ 

T 

T 

T 

1’ 

II 

’20072322242’ 

’1’ 

T 

T 

T 

T 

I 

’20072461614’ 

’0’ 

T 

T 

T’ 

T’ 

IV 

’20072491751’ 

’0’ 

’0’ 

T’ 

1’ 

’0’ 

IV 

’20072530149’ 

1’ 

T 

1’ 

1’ 

1’ 

II 
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'20072711338' 

'0' 

’1’ 

’1’ 

’1’ 

1’ 

IV 

'20073200313' 

1’ 

’1’ 

’1’ 

1’ 

1’ 

II 

'20073331900' 

1’ 

’1’ 

T 

T 

'O’ 

I/II 

'20073500809' 

'0' 

’1’ 

’0’ 

T 

’1’ 

II 

'20073530930' 

'1' 

’1’ 

’1’ 

'O’ 

’1’ 

IV 

'20080351701' 

'1' 

’1’ 

’0’ 

’0’ 

1’ 

II 

'20080390938' 

'1' 

’1’ 

'0' 

’1’ 

’0’ 

I 

'20080431250' 

'1' 

’0' 

’1’ 

’1’ 

’0’ 

II 

'20080451009' 

'0' 

’1’ 

’0’ 

'O’ 

’1’ 

I 

'20080520246' 

’1’ 

'0' 

’0’ 

’1’ 

'O’ 

I 
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3.  Magnetotelluric  imaging  of  the  Cascadia 
subduction  system  beneath  central  Washington 


Abstract 

In  this  chapter  we  present  images  from  the  inversion  of  magnetotelluric  data 
acquired  from  a  dense  line  of  60  broadband  and  20  long  period  instruments 
installed  in  central  Washington.  We  image  a  dipping  resistive  feature  in  the  west 
consistent  with  the  subducting  oceanic  slab.  We  also  image  two  conductors  that 
are  consistent  with  the  release  of  fluids  from  the  top  of  the  slab,  one  at  a  depth  of 
~  40  km  and  the  other  at  a  depth  of  80-90  km.  The  deeper  conductor  is  very 
conductive,  indicating  the  presence  of  substantial  fluid  and/or  melt,  and  rises  in  a 
column  to  another  conductive  feature  that  we  interpret  to  be  a  magma  chamber 
associated  with  Mt.  Rainier.  The  root  of  Mt.  Rainier  shows  up  as  a  resistive 
feature  that  extends  into  the  Western  Cascades,  giving  way  to  a  well  known 
conductive  feature  beneath  the  High  Cascades.  A  comparison  of  the  model  with 
the  seismic  image  from  Chapter  2  shows  that  the  rising  conductor  is  associated 
with  the  disruption  of  the  continental  Moho,  and  that  the  source  of  the  conductor 
at  80  km  depth  is  coincident  with  the  low  velocity  feature  in  the  migrated  seismic 
image  that  we  interpreted  as  fluids  or  melt.  While  the  image  is  consistent  with 
the  seismic  image  and  provides  substantial  new  information  on  its  own,  an 
assessment  of  the  model  itself  suggests  ways  in  which  it  might  be  improved,  a 
challenge  taken  up  in  the  final  chapter. 
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1.  Introduction 


Many  of  the  natural  hazards  associated  with  subduction  zones  are  intimately 
tied  to  the  release  of  fluids  bound  within  the  subducting  slab.  Intraslab 
earthquakes  have  been  linked  with  a  process  in  which  fluids  released  during  the 
transformation  of  hydrated  metabasalts  to  eclogite  in  the  upper  crust  increases 
pore  pressure,  thereby  reactivating  pre-existing  faults  (Kirby  et  al.,  1996,  Hacker 
et  al.,  2003b,  Preston  et  al.,  2003).  Slow  slip  and  low  frequency  tremor,  which 
have  been  tentatively  linked  with  the  periodicity  of  great  megathrust 
earthquakes  (Mazzotti  and  Adams,  2004),  can  be  explained  by  episodic  buildup 
and  release  of  fluid  pore  pressure  across  plate  boundaries  down-dip  of  the 
locked  zone  (Abers  et  al.,  2009,  Audet  et  al.,  2010).  Most  arc  magmas  are  also 
generally  believed  to  be  generated  by  flux  melting  that  occurs  when  dehydration 
reactions  in  the  oceanic  crust  and  uppermost  mantle  release  fluids  into  the 
continental  asthenosphere  (Ulmer  and  Trommsdorff,  1995,  Kirby  et  al.,  1996). 
Finally,  density  increases  associated  with  dehydration  have  been  linked  with 
changes  in  slab  dip  and  increased  stresses  within  the  slab  (Klemd  et  al.,  2011). 

Magnetotellurics  (MT)  is  an  attractive  method  to  apply  to  subduction  settings 
because  MT  is  particularly  sensitive  to  electrically  conductive  fluid  phases  and 
melt,  and  can  therefore  be  used  to  constrain  the  location  of  fluids  released  during 
dehydration  reactions.  Fluids  associated  with  dehydration  reaction  in  both  the 
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forearc  and  backarc  have  been  detected  using  MT  in  several  subduction  zones, 
including  Costa  Rica  (Worzewski  et  al.,  2011,  Brasse  et  al.,  2009)  and  Cascadia  to 
the  north  (Soyer  and  Unsworth,  2006)  and  south  (Wannamaker  et  al.,  1989,  Evans 
et  al.,  submitted)  of  the  present  study. 

Beneath  central  Washington,  the  Juan  de  Fuca  plate  is  subducting  at  a  rate  of 
approximately  35  mm/yr,  at  an  angle  20-35  degrees  north  of  east  (Fluck  et  al., 
1997,  Miller  et  al.,  2001).  The  dip  of  the  subducting  plate  is  approximately  10 
degrees  near  the  coast,  steepening  to  roughly  40  degrees  by  the  time  it  passes  the 
Cascade  Range  (McCrory  et  al.,  2006,  Roth  et  al.,  2008).  As  the  subducting  plate 
is  only  8-9  Mya  at  the  trench,  it  is  expected  to  release  much  of  its  fluid  budget  in 
the  forearc.  The  hydrated  metabasalt  of  the  upper  crust  will  eclogitize  at  a  depth 
of  approximately  40  km  in  Cascadia  (Hacker  et  al.,  2003b),  although  structural 
and  mineral  heterogeneities  will  spread  this  out  over  a  range  of  depths  (Kirby  et 
al.,  2006).  Hydrated  mantle  phases  such  as  serpentine  and  chlorite  will  remain 
stable  to  depths  of  more  than  90  km  in  Cascadia  (Hacker,  2003b)  and  have  been 
implicated  in  arc  magma  generation  (Hattori  and  Guillot,  2003). 

Although  consisting  of  fewer  than  30  high  peaks,  the  Cascade  Range  numbers 
more  than  2300  Quaternary  vents,  distributed  in  such  a  way  as  to  suggest 
segmentation  along-strike  (Hildreth,  2007).  The  most  abrupt  along-strike 
transition  occurs  in  central  Washington,  just  to  the  north  of  Mt.  Rainier  and  is 
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bisected  by  our  MT  station  profile  (Fig  1.)  To  the  north  of  our  line  there  is  a  gap 
of  120  km  to  the  Glacier  Peak  vent  cluster,  which  marks  the  southern  end  of  the 
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Fig  1:  Regional  map  showing  the  CAFE  MT  line  for  the  present  study  in  relation 
to  the  Cascadia  volcanic  arc.  Red  indicates  regions  where  forearc  and  arc  vents 
are  found,  and  back  arc  vent  regions  are  indicated  in  yellow.  P  is  the  Portland 
vent  field,  and  In  the  Indian  Heaven  vent  field.  (Adapted  from  Hildreth,  2007) 
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Garibaldi  Volcanic  Belt  (GVB).  The  GVB  consists  of  127  Quaternary  vents 
clustered  around  five  tall  stratovolcanoes  (Hildreth,  2007)  extending  in  a  line 
roughly  30  degrees  west  of  north.  In  contrast,  an  array  of  more  than  150 
Quaternary  vents  with  no  gaps  greater  than  10  km  between  them  extends 
roughly  10  degrees  west  of  south  for  170  km  to  Mt.  Hood  (Hildreth,  2007).  While 
most  of  the  volcanic  activity  takes  place  in  the  approximately  linear  axial 
volcanic  belt,  the  region  of  our  study  and  to  the  immediate  south  of  our  line  is 
notable  for  a  significant  presence  of  forearc  volcanism,  which  is  not  found 
elsewhere  in  Cascadia  (Hildreth,  2007).  This  includes  Mt.  Rainier  and  Mt.  St. 
Helens  in  addition  to  the  Indian  Heaven  and  Portland  vent  fields. 

The  CAFE  MT  experiment  consisted  of  a  dense  linear  array  of  60  wideband 
and  20  long-period  magnetotelluric  instruments  extending  approximately  280 
km  from  the  Washington  State  coast  north  of  Grays  Harbor,  passing  to  the  south 
of  Puget  Sound  and  Tacoma,  within  the  northern  boundary  of  Mt.  Rainier 
National  Park,  and  ending  about  25  km  east  of  Ellensburg  (Fig.2).  The  wideband 
stations  were  spaced  at  distances  of  5  km  or  less,  and  the  target  spacing  for  the 
long-period  stations  was  15  km.  The  CAFE  MT  experiment  was  designed  to  be 
collocated  with  the  earlier  CAFE  dense  array  of  passive  seismic  stations  (Fig.  2) 
in  order  to  utilize  constraints  from  migrated  images  obtained  during  that  phase 
(Abers  et  al.,  2009,  McGary  et  al.,  submitted). 
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Fig  2:  Map  showing  the  location  of  the  stations  for  the  CAFE  MT  and  seismic 
studies. 


2.  Fluids  in  subduction  settings 

The  composition  of  a  subducting  slab  can  be  idealized  into  four  distinct  layers. 
The  upper  crust  consists  of  hydrated,  fine-grained  pillow  basalts  and  sheeted 
dikes,  underlain  by  a  locally  hydrated,  coarse-grained,  low  porosity  gabbroic 
lower  crust.  Below  this  we  find  a  locally  hydrated  harzburgite/dunite  uppermost 
mantle,  and  below  that,  an  essentially  anhydrous  lower  mantle  layer  of  roughly 
lherzolite  composition.  The  extent  to  which  the  lower  crust  and  upper  mantle 
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are  locally  hydrated  can  vary  dramatically  depending  on  the  conditions  under 
which  they  were  formed  and  transported  to  the  subduction  zone.  This  is 
particularly  true  with  respect  to  faulting  history,  as  faults  can  provide  pathways 
for  fluids  to  infiltrate  these  layers. 

There  are  a  number  of  reactions  that  occur  within  a  subduction  system  that 
allow  for  the  release  of  fluids  from  the  descending  slab.  The  first  of  these  is  the 
conversion  of  hydrated  metabasalts  in  the  upper  crust  to  eclogite.  This  reaction 
is  largely  controlled  by  pressure  in  warm  subduction  settings  (Preston  et  al., 

2003,  Rondenay  et  al.,  2008),  and  occurs  at  depths  starting  at  -40-45  km  for 
Cascadia.  The  fine  grained  basalts  at  the  top  of  the  crust  eclogitize  first,  followed 
by  the  progressively  more  coarse  sheeted  gabbroic  dikes  below,  so  that  the 
reaction  actually  occurs  over  a  range  of  depths.  A  contribution  from 
eclogitization  of  locally  hydrated  lower  crust  is  possible  here  as  well,  extending 
the  depth  range  further. 

A  second  dehydration  reaction  that  occurs  within  the  subduction  setting 
involves  the  release  of  fluids  from  serpentinized  materials,  including  talc, 
brucite,  serpentine,  and  chlorite.  While  talc  and  brucite  are  not  stable  at  depth  in 
Cascadia,  the  dehydration  of  serpentine  is  controlled  largely  by  temperature,  and 
occurs  between  600-650°C  (Preston  et  al.,  2003,  Hattori  and  Guillot,  2003,  Audet 
et  al.,  2010),  while  chlorite  harzburgite  is  stable  to  temperatures  approaching 
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800°C  (Hacker  et  al.,  2003a).  Within  the  Cascadia  subduction  system,  this 
generally  corresponds  to  depths  of  80-95  km  (Hacker  et  al,  2003). 

The  introduction  of  serpentinized  material  to  such  depths  in  the  subduction 
system  can  occur  in  a  variety  of  ways.  The  antigorite  serpentine  and  chlorite  that 
form  at  the  bottom  of  the  mantle  wedge  can  be  dragged  to  depth  by  the 
movement  of  the  subducting  slab  (Hattori  and  Guillot,  2003,  Kawakatsu  and 
Watada.,  2007),  but  the  thermal  profile  for  Cascadia  excludes  this  possibility 
(van  Keken  et  al.,  2011).  Serpentinized  mantle  rocks  have  also  been  identified  in 
oceanic  lithosphere  at  slow  and  intermediate  rate  spreading  ridges,  but  only  in 
regions  where  the  axial  valleys  are  well  developed  (Cannat,  1993),  a  condition 
that  is  not  met  for  the  southern  Juan  de  Fuca  ridge  (Canales  et  al.,  2005,  Carbotte 
et  al.,  2006). 

For  Cascadia,  any  hydration  of  mantle  rocks  is  likely  the  result  of  fluid 
intrusion  that  occurs  due  to  the  reactivation  of  the  normal  faults  formed  at  the 
spreading  center.  These  faults  have  been  shown  to  reactivate  and  provide  fluid 
transport  pathways  into  the  lower  crust  and  uppermost  mantle  during  bending 
at  the  trench  (Ranero  et  al.,  2005,  Key  et  al.,  2012),  but  for  Cascadia,  the  shallow 
dip  angle  (McCrory  et  al.,  2006)  and  the  hydraulic  barrier  presented  by  the  thick 
sediment  layer  at  the  trench  (Divens  et  al.,  retrieved,  2012)  are  formidable 
obstacles  to  hydration  in  this  manner.  Hydrated  mantle  rocks  have  been 
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identified  in  Cascadia  as  a  consequence  of  reactivation  of  faults  within  potential 
weak  zones  corresponding  to  wake  propagation  (Nedimovic  et  al,  2009),  and 
have  been  implicated  as  the  cause  of  the  2001  magnitude  6.8  Nisqually 
earthquake  (Kao  et  al.,  2008). 

The  coarse  grained,  non-porous  gabbros  of  the  lower  crust  may  also  transform 
directly  to  eclogite,  a  reaction  which  can  occur  at  depths  of  85  to  110  km  for 
Cascadia  (John  and  Schenk,  2003),  but  this  reaction  tends  to  be  very  sluggish  in 
the  absence  of  free  water  (Hacker  et  al.,  2003b).  A  key  catalyst  is  fluid  access  to 
pathways  within  the  lower  crust,  therefore  it  is  expected  that  dehydration 
reactions  of  serpentine  and  chlorite  in  the  underlying  upper  mantle  would 
accelerate  this  process  (John  and  Schenk,  2003). 

Fluids  released  during  these  hydration  reactions  migrate  upwards  into  the 
overlying  mantle  wedge,  where  their  presence  lowers  the  melting  temperature  of 
the  mantle  rock  potentially  leading  to  flux  melting.  The  combination  of  the 
thermal  conditions  and  the  presence  of  fluids  are  thought  to  control  the 
formation  and  distribution  of  arc  volcanism  (Grove  et  al.,  2009). 

3.  Data  and  Processing 
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Time  series  data  from  20  long-period  and  60  wideband  MT  instruments  were 
collected  for  the  CAFE  MT  experiment.  The  long-period  data  were  collected 
using  NIMS  instruments  with  each  station  in  place  for  typically  three  weeks.  The 
wideband  data  were  collected  by  Quantec  Geoscience  Inc.  using  Phoenix 
Geosystems  instrumentation,  with  a  typical  recording  interval  of  15  hours  for 
each  site.  The  280  km  long  array  through  central  Washington  was  designed  to  be 
roughly  collocated  with  the  earlier  CAFE  passive  seismic  experiment. 

The  time  series  data  from  each  long-period  station  was  visually  inspected  for 
breaks,  trends,  and  signal  to  noise  ratio,  then  windowed  with  data  from  two 
separate  magnetic  remote  reference  sites.  One  remote  reference  for  each  station 
was  selected  from  distant  CAFE  MT  stations,  and  the  second  concatenated  from 
Earthscope/US  Array  stations  in  Nevada  and  California  that  were  operated 
concomitantly.  The  minimum  distance  between  the  primary  station  and  the 
second  remote  was  greater  than  500  km,  giving  confidence  that  the  noise 
between  the  stations  was  not  likely  to  be  correlated  to  any  significant  degree. 
Dual  remote  references  were  also  used  for  the  wideband  data,  with  one  remote 
reference  located  approximately  30  km  east  of  Mt.  Rainier,  and  a  distant  remote 
reference  located  at  Buena  Vista  Valley,  Nevada. 

The  long-period  time  series  data  were  processed  into  impedance  tensors  using 
the  robust  bounded  influence  BIRRP  algorithm  (Chave  and  Thompson,  2003, 
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2004).  The  wideband  data  were  processed  following  the  procedure  of  Egbert  and 
Booker  (1996).  The  long-period  transfer  functions  provided  useful  data  from  10 
seconds  up  to  several  thousand  seconds,  and  the  wideband  from  less  than  one 
second  up  to  several  hundred  seconds.  The  data  were  collected  during  a  low 
solar  cycle,  and  were  also  affected  by  local  noise,  particularly  near  the  cities  of 
Olympia  and  Tacoma,  so  the  period  range  of  the  useful  data  was  not  as  long  as  it 
might  have  been  otherwise.  Samples  of  both  the  long-period  and  wideband 
transfer  functions  are  shown  below  (Fig.  3).  For  the  inversions,  these  data  were 
supplemented  with  transfer  functions  from  four  stations  that  were  part  of  the 
Earthscope  US  Array  in  Washington. 


CAFE  LP  09  CAFE  WB  37 


Fig  3:  Sample  transfer  functions  for  a  long  period  station  (FP  09)  and  a  wide  band 
station  (WB07) 
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After  correcting  the  data  for  declination,  dimensionality  and  regional  strike 
direction  were  evaluated  for  the  data  in  a  variety  of  ways,  including  phase  tensor 
analysis  (Caldwell  et  al.,  2004),  Bahr  skew  analysis  (Simpson  and  Bahr,  2005), 
and  multisite,  multifrequency  Groom-Bailey  tensor  decomposition  using  the 
Strike  software  package  (McNeice  and  Jones,  2001,  Groom  and  Bailey,  1989).  The 
phase  tensor  ellipses  (Fig.  4)  show  a  consistent  pattern  for  periods  between  20 
and  2000  seconds,  with  a  near  north-south  strike  for  stations  west  of  Mt.  Rainier, 
shifting  to  a  slightly  clockwise  strike  to  the  east  of  Mt.  Rainier,  with  a  more 
pronounced  clockwise  change  for  stations  towards  the  eastern  end  of  the  profile. 
This  is  consistent  with  previous  strike  analysis  conducted  on  a  set  of  Earthscope 
LP  stations  along  a  line  at  roughly  46.5  N  (Patro  and  Egbert,  2008). 

Bahr  (phase-sensitive)  skew  analysis  provides  a  justification  for  a  2-D  analysis 
of  the  data  between  the  same  periods,  with  exceptions  beneath  stations  LP-32 
and  LP-36  (immediately  west  of  Mt.  Rainier)  for  periods  greater  than  400 
seconds,  and  beneath  stations  LP-44  (just  east  of  Mt.  Rainier)  and  LP-54  (in  a 
river  valley  just  north  of  Ellensburg)  for  periods  less  than  150  seconds.  It  has 
been  shown  that  while  we  can  consider  these  exceptions  as  exhibiting  3-D 
characteristics,  we  cannot  assume  that  the  rest  of  the  data  can  be  considered  2-D 
based  strictly  on  phase-sensitive  skew  analysis  (Ledo  et  al.,  2002). 
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The  Groom-Bailey  decomposition  technique  separates  from  the  impedance 
tensor  determinable  matrices  representing  twist,  shear,  and  anisotropy,  leaving 
only  distortion  elements  that  affect  both  the  amplitude  and  phase  of  the  electric 


Fig  4:  A  map  of  the  Pacific  Northwest  showing  the  CAFE  long  period  stations  as 
blue  circles.  The  ellipses  show  the  average  direction  of  the  phase  tensor-ellipse 
for  periods  between  100  and  1340  seconds.  (Caldwell  et  al.,  2004). 


field,  and  performs  statistical  analysis  of  the  validity  of  the  assumption  of 
galvanic  distortion  (Groom  and  Bailey,  1989,  Jones,  2012).  The  Strike  algorithm 
(McNeice  and  Jones,  2001)  is  an  extension  to  this  decomposition  that  allows 
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determination  of  a  best  fitting  strike  direction  for  multiple  sites  and  over 
multiple  frequencies.  For  the  CAFE  long-period  data,  a  regional  strike  direction 
of  seven  degrees  west  of  north  provided  the  best  chi-squared  fit  for  the  long- 
period  data  for  periods  of  100-1350  seconds,  within  confidence  intervals  when 
stations  LP-32  and  LP-36  were  excluded.  A  strike  direction  of  due  north  also 
produced  an  acceptable  fit,  and  was  used  in  the  inversions  (i.e.  the  data  were 
decomposed  but  not  rotated  using  the  Strike  algorithm).  The  strike  angles  for  the 
entire  set  of  decomposed  transfer  functions  are  displayed  in  the  form  of  a  rose 


diagram  (Fig.  5). 
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Fig  5:  Rose  diagram  showing  overall  strike  directions  with  the  color  code 
reflecting  the  degree  of  anisotropy  as  determined  using  the  STRIKE  algorithm  for 
the  CAFE  data  set  (McNeice  and  Jones,  2001). 
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4.  Modeling 


Our  2-D  inversion  model  assumes  a  primary  strike  direction  of  due  north, 
consistent  with  the  decomposition  efforts  described  in  the  previous  section. 

To  generate  our  regularized  2-D  models  of  resistivity  structure,  we  used  the  non¬ 
linear  conjugate  gradient  inversion  algorithm  WinGLink  (Rodi  and  Mackie, 

2001).  The  inversions  were  run  starting  with  a  half-space  of  100  Qm,  except  for 
the  ocean,  which  was  fixed  at  0.33  Qm  with  an  extent  defined  by  local 
bathymetry.  The  inversion  mesh  consisted  of  116  rows  and  234  columns.  The 
row  height  was  very  fine  near  the  surface,  increasing  gradually  with  depth, 
while  the  column  width  was  maintained  to  be  generally  uniform  as  allowed  by 
the  station  spacing. 

The  inversion  process  is  iterative  and  attempts  to  minimize  a  Tikhonov 
objective  function,  which  can  be  expressed  as: 

'P(m)  =  (d  -  F(m))T  V1  (d-F(m))+A  mT  LT  Lm  (1) 

where  F(m)  is  the  forward  modeling  function  and  d  the  data  vector  (with  each  d; 
representing  the  log  apparent  resistivity  or  phase  for  a  given  polarization, 
observations  site,  and  frequency),  V  is  a  matrix  that  represents  the  variance  of  the 
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error  vector,  and  A  is  the  (generally  global)  regularization  parameter  that 
balances  the  misfit  term  with  the  smoothness  term  LTLm,  which  for  our 
inversions  is  defined  by  the  integral  of  the  square  of  the  Laplacian.  In  the 
WinGLink  implementation  of  the  inversion  code,  the  tradeoff  between 
smoothness  and  misfit  is  defined  by  the  value  x  (instead  of  A),  and  the  optimal 
value  is  determined  by  plotting  rms  misfit  against  model  roughness.  For  our 
inversions  a  x  value  of  3.3  was  deemed  to  be  optimal.  Two  other  parameters,  a 
and  [3,  define  the  relationship  of  the  smoothness  parameter  in  the  vertical  as 
compared  to  the  horizontal,  and  the  way  in  which  the  smoothness  changes  with 
depth  respectively.  A  range  of  values  from  0.8  to  1.8  for  a  and  1.0  to  2.0  for  [3 
were  evaluated,  with  final  selected  values  of  a=1.5  and  (3=1.7.  While  small 
changes  in  a  and  (3  have  been  shown  to  produce  striking  differences  in  structure 
in  some  cases  (Matsuno  et  al.,  2011),  the  changes  that  we  saw  in  structure  while 
varying  these  parameters  was  quite  minimal. 

Transfer  functions  for  a  2-D  MT  inversion  consist  of  impedance  and  phase 
information  over  a  range  of  frequencies,  for  two  perpendicular  modes,  the 
transverse  magnetic  (TM)  and  transverse  electric  (TE).  We  can  additionally 
utilize  the  tipper  function,  which  relates  the  horizontal  magnetic  field 
information  to  the  vertical  magnetic  field.  The  relationship  between  these 
parameters  can  be  expressed  as  in  equation  1, 
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where  Ex  and  Ey  are  the  two 

perpendicular  horizontal  components  of  the  electric  field,  and  Bx,  By,  and  Bz  are 
the  two  perpendicular  horizontal  components  and  the  vertical  component  of  the 
magnetic  field.  The  Z  components  represent  the  complex  impedences  (apparent 
resistivity  and  phase)  with  the  caveat  that  Zxy  and  Zyx  are  zero  with  the  correct 
rotation  of  x  and  y  for  an  ideal  two  dimensional  subsurface.  The  components  of 
the  tipper  function  are  shown  as  Tzx  and  Tzy. 

We  generated  models  that  inverted  both  the  TM  and  TE  modes,  with  error 
floors  of  5%  and  15%  respectively,  and  a  set  of  models  using  only  the  TM  mode, 
which  has  been  shown  to  be  less  influenced  by,  and  therefore  more  robust  to  the 
effects  of  3D  structure  (Wannamaker  et  al.,  1989).  The  tipper  function  was 
included  in  both  sets  of  inversions. 

We  also  evaluated  the  extent  to  which  the  primary  structures  in  our  models 
were  supported  and/or  required  by  the  data.  This  was  done  in  a  number  of 
ways,  such  as  assessing  how  the  structures  were  affected  by  parameter 
modifications  such  as  changes  to  the  variables  controlling  global  smoothness  (i.e. 
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t,  a,  and  [3  described  above),  and  comparing  the  resulting  models  to  the  data 
pseudo  sections  to  give  us  a  better  understanding  of  how  data  from  each  station 
might  be  affecting  the  inversion. 

Additionally,  manual  editing  techniques  were  utilized.  This  involved 
removing  or  altering  a  feature  in  the  resultant  model  and  then  allowing  the 
inversion  process  to  seek  a  solution  optimally  close  to  this  altered  model.  If  the 
feature  was  restored  by  the  inversion,  it  was  taken  to  be  required  by  the  data.  If 
the  inversion  was  able  to  achieve  a  misfit  comparable  to  the  original  misfit 
without  restoring  the  feature,  the  structure  in  question  was  determined  to  be 
allowed  but  not  required  by  the  data. 

Finally,  in  order  to  assess  the  inverted  image  as  a  whole,  we  generated  a  set  of 
synthetic  station  data  from  the  inverted  images,  and  then  performed  a  series  of 
inversions  using  the  synthetic  data  and  starting  from  our  100  Cm  half-space 
model  with  fixed  resistivity  for  water.  We  should  expect  that  if  the  inverted 
image  were  completely  correct,  the  inversion  of  the  synthetic  data  would 
faithfully  reproduce  the  image,  whereas  variations  between  the  images  can 
provide  important  information  regarding  the  limitations  of  the  inversion.  This 
idea  will  be  discussed  in  more  detail  in  section  6. 

5.  Results 
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Figures  6a  (top)  and  6b  (bottom)  show  our  primary  inversion  models.  The 
image  in  figure  6a  was  generated  using  only  the  TM  mode  along  with  the  tipper, 
while  the  image  in  figure  6b  also  included  the  TE  mode  in  the  inversion.  The 


40 

|. 


80 

120 

50  100  150  200  250 

Distance  from  coast  (km) 


Figure  6:  This  figure  shows  the  primary  inversion  images  for  the  CAFE  data.  The 
top  image  was  generated  using  just  the  TM  mode  and  the  tipper  function  data, 
while  the  bottom  image  included  the  TE  mode  data 


two  inversions  produce  images  that  are  very  similar,  a  conclusion  supported  by 
generating  pseudo-sections  from  each  of  the  inversion  models,  and  differencing 
them  to  see  the  contrast.  Differences  in  the  apparent  resistivities  in  the  models 
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never  exceeded  10%,  and  differences  in  the  phases  remained  less  than  2° 
everywhere.  Despite  this  similarity,  the  rms  misfit  for  the  TM/TE/tipper  model  is 
3.03,  significantly  higher  than  that  for  the  TM/tipper  model  at  1.96.  This  suggests 
that  both  inverted  images  were  driven  primarily  by  the  TM  and  tipper  data. 


Figure  7  shows  plots  of  rms  misfit  value  by  station  for  each  of  the  inversion 
models.  For  the  most  part,  the  shapes  of  the  plots  are  similar,  with  the  rms  misfit 
values  staying  fairly  close  between  models  (though  they  are  consistently  higher 
for  the  TM/TE/tipper  inversion).  The  clear  exceptions  to  this  are  between 
stations  6-10,  corresponding  to  the  western  part  of  the  Willapa  Valley,  and 
between  stations  36  and  42,  corresponding  to  the  vicinity  of  Mt.  Rainier,  where 
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the  increasing  misfit  of  the  TM/TE/tipper  inversion  is  contrasted  by  a  steady  or 
decreasing  misfit  for  the  TM/tipper  only  inversion.  While  it  might  be  tempting 
to  invoke  3D  effects,  this  is  not  consistent  with  the  Bahr  skew  and  strike  analyses. 
This  discrepancy  has  more  to  do  with  our  standard  choice  of  a  globally  smooth 
half  space  model,  which  will  be  examined  in  detail  in  Chapter  4.  For  the  time 
being  we  will  set  this  aside,  and  continue  with  our  interpretation  of  the 
TM/TE/tipper  inversion  model. 

Shallow  features 

There  are  three  shallow  conductors  in  the  inversion  cross-sections  (Fig.  8)  that 
correspond  to  sediment  basins  along  the  transect.  The  first  (Al)  occurs  in  the 
Willapa  Valley,  between  the  coastal  Grays  Harbor  region  and  the  uplifted  Black 
Hills.  The  second  (A2)  corresponds  with  the  quaternary  glacial  deposits 
(Schuster,  2007)  in  the  forearc  basin,  and  the  third  (A3)  with  the  quaternary 
sediments  in  the  backarc  basin.  These  shallow  conductors  are  typical  of  those 
previously  found  during  MT  work  in  other  subduction  settings  (Soyer  and 
Uns worth,  2006,  Brasse  et  al.,  2009,  Evans  et  al.,  submitted). 

A  highly  resistive  crust  (B)  measuring  several  hundreds  to  a  few  thousand 
Qm  underlies  Rainier  National  Park,  consistent  with  the  presence  of  the  dry 
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pyroxene-andesite-dacite  complex  that  makes  up  most  of  the  volcanic  basement 
(Simpson  and  Bahr,  2005,  Hildreth,  2007).  Resistivities  in  the  Western  Cascades 


Figure  8:  The  upper  60  km  of  the  conductivity  structure  for  the  CAFE  line.  The 
red  triangle  marks  the  location  of  Mt.  Rainier  just  to  the  south  of  the  profile.  Al- 
A3  mark  conductors  associated  with  sedimentary  basins.  B,  Cl,  and  C2  are 
resistive  features  associated  with  Mt.  Rainier  and  the  Western  Cascades.  D  is  a 
conductive  feature  beneath  the  High  Cascades.  The  resistor  E  is  associated  with 
the  crustal  forearc,  and  we  interpret  the  strong  conductor  G  as  a  magma  chamber 
feeding  Mt.  Rainier. 


(Cl)  are  also  in  the  thousands  of  Cm ,  but  drop  off  quickly  into  the  High 
Cascades  (D),  where  they  range  from  roughly  10-100  Qm.  This  transition  in 
resistivity  from  the  Western  Cascades  to  the  High  Cascades  has  been  shown  to 
extend  much  of  the  length  of  the  Cascade  Range  (Jiracek  and  Curtis,  1989, 
Wannamaker  et  al.,  1989,  Patro  and  Egbert,  2008),  and  is  thought  to  indicate  the 
presence  of  interconnected  fluids  (Wannamaker  et  al.,  1989).  This  conductive 
feature  dips  towards  the  east,  and  appears  to  extend  potentially  all  of  the  way  to 
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the  upper  mantle,  an  observation  that  has  also  been  previously  reported  (Patro 
and  Egbert,  2008).  Referred  to  as  the  Southwest  Washington  Cascades 
Conductor  (SWCC),  it  is  described  as  an  upper  crustal  conductive  anomaly  just 
west  of  the  Cascades,  and  it  has  been  suggested  that  the  deeper  parts  of  this 
conductor  may  have  a  distinct  cause  such  as  fluids  associated  with  subduction 
and  arc  magmatism  (Patro  and  Egbert,  2008). 

The  crust  in  the  forearc  region  (E)  is  generally  quite  resistive,  from  several 
hundred  to  more  than  a  thousand  Cm,  not  quite  as  resistive  as  might  be  expected 
for  a  tectonically  inactive  upper  crust  (Nesbitt,  1993,  Gough,  1986),  but  consistent 
with  what  has  been  observed  in  other  subduction  settings  (Brasse  et  al.,  2009).  In 
Costa  Rica,  this  discrepancy  was  attributed  to  fluid  input  derived  from 
dehydration  of  the  downgoing  slab,  and  this  conclusion  is  very  strongly 
supported  in  the  current  case  as  well.  We  will  examine  this  in  more  detail  when 
we  consider  the  inversion  image  within  the  context  of  previous  passive  seismic 
work  in  section  5.2. 

The  most  striking  feature  within  the  crust  is  the  strong  conductor  (G)  located 
at  depths  from  around  10  km  roughly  below  Mt.  Rainier.  While  much  of  the 
volcanic  rock  that  makes  up  Mt.  Rainier  is  pyroxene-andecite-dacite,  inclusions 
of  more  basaltic  andesite  are  common  and  most  likely  represents  the  main  parent 
magma  in  the  crustal  reservoir  (Hildreth,  2007).  These  samples  typically  give  an 
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Si02  content  of  -55% wt  and  Na20  content  of  ~3.5%wt  (Stockstill  et  al.,  2002). 


Assuming  a  typical  water  contribution  of  ~4%wt  and  a  temperature  of  1275  K 
consistent  with  the  solidus  for  this  magma,  we  would  expect  a  resistivity  of  -0.7 
Qm  for  a  pure  melt  at  15  km  depth  (Pommier  and  LeTrong,  2007).  The  lowest 
resistivities  in  the  model  are  -0.85  Qm,  suggesting  substantial  pockets  of  high 
melt  fractions  (close  to  100%).  The  specific  distribution  of  conductivities  within 
the  anomaly  is  consistent  with  the  existence  of  two  pockets,  a  primary  pocket  to 
the  west  of  Mt.  Rainier  at  a  depth  of  15  km,  and  a  secondary  source  at  10  km 
depth  more  directly  beneath  the  volcano.  The  mixed  and  zoned  eruptive  history 
of  Mt.  Rainier  has  led  others  to  argue  for  such  multi-level  reservoirs  (Venezky 
and  Rutherford,  1997,  Stockstill  et  al.,  2002).  The  primary  reservoir  clearly 
appears  to  be  connected  to  a  source  in  the  mantle;  we  will  consider  this  in  more 
detail  in  section  5.2. 

To  the  east  of  this  conductor,  the  resistor  underlying  the  Western  Cascades 
can  be  seen  to  extend  to  or  be  connected  to  another  resistive  feature  (C2)  that 
extends  to  depths  of  25  km.  The  resistivity  of  this  feature  is  on  the  order  of  a  few 
thousand  Qm,  consistent  with  either  tectonically  inactive  upper  crust  (Nesbitt, 
1993,  Gough,  1986)  or  intruded,  cool,  mafic  igneous  rock  (Simpson  and  Bahr, 
2005). 
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Deeper  features 


Previous  work  has  shown  that  MT  data  is  not  always  particularly  sensitive  to 
the  structure  within  a  subducting  slab,  sometimes  showing  a  varying  resistive 
feature  that  roughly  coincides  with  the  top  of  the  slab  (Soyer  and  Unsworth, 

2005,  Evans  et  al.,  submitted),  and  sometimes  returning  no  commensurate 
resistive  region  at  all  (Matsuno  et  al.,  2011).  The  unconstrained  model  in  this 
study  is  an  intermediate  case,  where  a  resistive  feature  is  present  (H  in  figure  9b), 
but  does  not  strictly  coincide  with  the  top  of  the  slab,  and  is  also  thicker  than 
would  be  expected  based  on  simple  thermal  considerations.  This  result  is  similar 
to  those  seen  previously  in  Central  America  (Brasse  et  al.,  2009),  and  the 
resistivity  of  more  than  1000  Cm  is  consistent  with  that  expected  for  the  slab 
(Bedrosian,  2007).  We  will  explore  some  of  the  variables  that  affect  the 
limitations  of  MT  as  regards  defining  the  top  of  the  resistive  slab  in  Chapter  4,  as 
well  as  some  methods  for  minimizing  these  limitations,  but  for  the  time  being  we 
will  use  the  results  from  the  collocated  passive  seismic  experiment  (see  Fig.  2  for 
stations)  in  Chapter  2  to  constrain  the  location  of  the  descending  oceanic  slab. 

Figure  9a  (top)  shows  the  image  generated  using  GRT  migration  of  passive 
seismic  data,  as  well  as  the  interpretation  of  that  image  discussed  in  Chapter  2. 
Figure  9b  (bottom)  shows  the  TM/TE/tipper  magnetotelluric  inversion  image 


91 


from  Fig.  5b,  but  with  the  seismic  interpretation  superimposed  to  provide  a 
framework  for  the  discussion  of  the  processes  illuminated  by  the  distribution  of 
conductivities. 

Starting  at  a  depth  of  approximately  40  km  and  intersecting  roughly  with  the 
top  of  the  descending  slab,  we  find  the  corner  of  a  v-shaped  conductor  (F)  with  a 
resistivity  consistently  around  80  Qm.  This  is  similar  to,  though  slightly  more 
resistive  than,  features  found  at  similar  depths  in  other  subduction  zone  studies 
(Brasse  et  al.,  2011,  Evans  et  al.,  submitted),  where  it  has  been  interpreted  as 
migrating  fluids  that  have  been  released  from  the  subducting  oceanic  crust.  This 
interpretation  is  also  strongly  supported  by  the  seismic  work,  as  a  velocity 
increase  in  the  crust  suggests  a  fluid  release  due  to  the  transition  of  hydrated 
metabasalts  to  eclogite  at  this  location. 

The  80  Qm  resistivity  is  consistent  with  a  0.4%wt  of  saline  fluids  (consistent  with 
seawater)  given  good  connectivity  (Pommier  and  LeTrong,  2007),  with  a  tradeoff 
between  higher  fluid  content  and  poorer  connectivity  also  yielding  acceptable 
results.  The  westward  extension  of  the  conductor  follows  the  top  of  the  slab  for 
~20  km  before  rising  directly  upward,  apparently  connecting  with  the  eastern 
edge  of  the  Willapa  Bay  sediments  conductor.  The  eastern  extension  rises 
approximately  20  km  over  a  horizontal  distance  of  40  km,  and  possibly 
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continuing  along  that  path  to  eventually  connect  with  the  conductor  (G)  that  we 


have  interpreted  as  the  magma  chamber  below  Mt.  Rainier. 
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Figure  9:  Figure  9a  (top)  shows  the  migrated  passive  seismic  image  from 
Chapter  2,  including  the  interpretation.  The  MT  model  is  shown  in  figure  9b 
(bottom).  The  resistors  H  and  I  are  associated  with  the  subducting  slab  and  the 
cold  nose  of  the  mantle  wedge,  respectively.  The  conductive  feature  F  shows  the 
release  of  fluids  from  eclogitization  of  the  upper  crust,  and  the  strong  conductor  J 
the  fluids  most  likely  with  a  component  of  flux  melting  associated  with  the  deep 
fluid  release. 


Following  the  slab  crust  down  to  a  depth  of  ~50  km,  we  encounter  a  feature  (I) 
with  a  resistivity  between  300-500  Cm.  Similar  features  have  been  encountered 
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in  earlier  subduction  zone  MT  studies  (Soyer  and  Unsworth,  2005,  Brasse  et  al., 
2009).  The  conductivities  are  consistent  with  either  mantle  minerals  or  antigorite 
serpentine  with  a  well  connected  fluid  phase  of  less  that  l%wt  (Reynard  et  al., 
2001,  Pommier  and  LeTrong,  2007),  but  also  with  dry  peridotite  at  the 
temperatures  expected  in  the  wedge.  This  is  consistent  with  the  seismic  image, 
as  this  resistor  corresponds  to  the  section  of  mantle  wedge  that  does  not  get 
significant  fluid  input  from  slab  dehydration. 

Following  the  top  of  the  slab  to  a  depth  of  ~80  km,  we  encounter  a  large  and 
rather  intense  conductor  (J),  with  resistivities  as  low  as  2  Qm.  The  bottom  of  this 
conductor  appears  to  coincide  with  the  top  of  the  slab  and  connect  with  the 
shallow  conductor  (G)  that  we  interpreted  as  a  magma  chamber.  The  implication 
is  that  this  conductor  represents  melt  caused  by  flux  melting  as  a  consequence  of 
fluid  release  from  a  depth  of  80-100  km,  and  that  this  melt  provides  the  source 
for  the  magma  chamber  that  supplies  Mt.  Rainier.  A  basaltic  andesite  magma  of 
the  type  inferred  to  be  the  source  magma  for  Mt.  Rainier  (Stockstill  et  al.,  2002, 
Hildreth,  2007),  at  a  temperature  of  1300  K  and  pressure  of  2.4  MPa  would  have  a 
resistivity  as  low  as  0.6  Qm  for  a  pure  melt  with  a  4%wt  water  content  (Pommier 
and  LeTrong,  2007),  which  is  consistent  with  the  most  conductive  area  of  the 
feature,  and  suggests  partial  melting  in  the  upward  extension  that  feeds  the 
crustal  magma  chamber.  The  most  conductive  zones  of  the  anomaly  are  also 
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coincident  with  the  low  velocity  region  in  the  seismic  image  that  we  interpreted 
as  a  fluid  melt  phase  based  on  independent  seismic,  hypocenter,  and  petrological 
evidence  in  Chapter  2.  The  conductive  feature  appears  to  be  entirely  consistent 
with  the  interpretation  in  that  chapter,  specifically  that  the  dehydration  of 
serpentine  and  chlorite  harzburgite  in  the  upper  mantle  releases  substantial 
fluid,  triggering  a  chain  reaction  which  transforms  the  metastable  gabbros  in  the 
lower  crust  into  eclogite.  The  released  fluids  migrate  into  the  overlying 
continental  mantle,  leading  to  flux  melting  and  the  formation  of  a  significant 
magma  source.  This  magma  source  feeds  the  crustal  level  chamber  or  chambers 
that  in  turn  provide  source  magma  for  Mt.  Rainier.  The  evidence  from  the  three 
seismic  lines  and  the  hypocenter  locations  of  deep  earthquakes  (addressed  in 
detail  in  Chapter  2)  suggest  that  this  mechanism  is  not  ubiquitous  along  the 
strike  of  the  Cascades.  This  is  not  unexpected,  however,  as  Mt.  Rainier  is  a 
forearc  volcano  some  50  km  to  the  west  of  the  active  Cascade  arc,  meaning  that 
our  interpretation  would  be  supported  should  similar  features  be  discovered  in 
the  vicinity  of  Mt.  St.  Helens,  the  only  other  major  forearc  volcano  in  the 
Cascades. 
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6.  Assessment  of  the  model 


Each  of  the  key  features  discussed  in  section  5  were  subject  to  a  number  of 
sensitivity  tests  in  order  to  determine  if,  within  the  context  of  the  remaining 
structures,  they  were  required  by  the  data.  This  was  done  by  removing  or 
changing  a  feature  and  then  allowing  the  inversion  code  to  seek  a  solution 
optimally  close  to  this  altered  model.  If  the  solution  was  unable  to  achieve  a 
similar  or  better  misfit  or  if  the  feature  was  simply  restored  by  the  inversion,  we 
concluded  that  the  original  feature  is  required  by  the  data. 

We  also  generated  a  set  of  synthetic  data  from  our  final  model,  and  then 
performed  a  series  of  inversions  using  the  synthetic  data  in  place  of  the  real 
processed  data.  The  inversions  were  produced  using  the  same  basic  100  Cm 
halfspace  with  a  0.33  Cm  ocean  determined  by  bathymetry  that  we  used  for  the 
original  models.  If  our  original  model  were  a  completely  correct  representation 
of  the  subsurface,  we  should  expect  that  the  inverted  image  using  the  synthetic 
data  would  faithfully  reproduce  the  original  model.  Of  course,  the  model  will 
never  be  perfect,  and  in  this  case  we  had  reason  to  suspect  places  where  the 
reproduction  would  highlight  flaws  given  the  mismatch  between  the  rms  misfits 
for  the  TM/tipper  and  TM/TE/tipper  inversions. 
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The  model  generated  from  the  synthetic  data  is  compared  to  the  original 
inversion  model  in  figure  10,  and  at  a  glance  it  is  clear  that  all  of  the  primary 
features  of  the  original  model  are  reproduced.  This  gives  us  confidence  that  our 
model  is,  at  least  in  general,  largely  correct.  Despite  this,  there  are  some  clear 
differences  between  this  model  and  our  primary  model. 

First,  the  conductor  that  we  interpreted  as  the  magma  chamber  beneath  Mt. 
Rainier  is  much  less  defined  in  the  synthetic  reproduction.  The  implication  is 
that  details  in  the  structure  of  a  conductive  feature  can  be  lost  during  the 
inversion  process.  This  suggests  that  the  actual  distribution  of  the  conductor  in 
the  subsurface  might  be  more  complicated  than  that  in  our  model,  an  idea  that 
we  will  explore  in  more  detail  in  Chapter  4. 

The  conductor  below  the  High  Cascades,  while  still  present  in  the  synthetic 
model,  shows  not  only  a  loss  of  detail,  but  the  evidence  for  the  connection 
between  this  conductor  and  the  conductive  mantle  to  the  east  is  much  less 
convincing. 

The  V-shaped  conductor  in  the  west  that  we  interpreted  as  fluids  released 
during  the  eclogitization  of  the  hydrated  basalts  in  the  upper  crust  has  also  lost 
its  shape,  smearing  out  along  a  horizontal  line  at  a  depth  of  ~50  km.  The  shape 
of  the  underlying  resistive  feature  has  also  been  altered.  We  would  expect  that 
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the  actual  conductivity  structure  of  fluids  being  released  from  the  upper  crust  of 
a  slab  to  be  represented  by  a  conductor  rising  from  the  top  of  a  strong  resistive 
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Figure  10:  This  figure  compares  the  original  inversion  model  (top)  with  the 
model  using  synthetic  data  (bottom).  The  synthetic  data  is  generated  using  the 
original  inversion  model  as  the  forward  model. 


feature,  and  while  the  inversion  is  consistent  with  this,  the  location  of  the 
interface  is  not  clearly  identified  in  the  inversion. 

Finally,  the  large  conductor  that  we  interpret  as  the  fluid  release  from 
serpentinization  at  ~  80-100  km  depth  is  both  more  compact  and  less  intense  in 
the  model  generated  from  synthetic  data.  It  is  also  worth  noting  that  in  this 
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model,  the  resistor  to  the  west  that  would  represent  the  slab  material  has 
retreated  from  the  conductor.  This  is  another  situation  where  we  would  expect  a 
strong  conductor  and  a  resistor  to  be  essentially  in  contact,  and  again  the 
inversion  does  not  clearly  identify  the  location  of  the  expected  interface.  Instead, 
the  smoothness  parameter  forces  the  resistor  and  conductor  to  separate,  and  in 
releasing  the  contrast,  weakens  the  intensity  of  both  features.  We  also  notice  that 
the  most  intense  region  of  the  conductor  in  both  inversions  is  centrally  located 
within  the  conductor,  when  we  would  expect  to  find  it  closer  to  the  top  of  the 
subducting  slab. 

In  Chapter  4,  we  will  examine  the  mechanisms  underlying  each  of  these 
phenomena,  and  also  look  for  ways  to  apply  this  information  within  the 
inversion  process  to  more  faithfully  reproduce  the  underlying  conductivity 
structure.  In  the  end,  we  will  generate  an  inversion  model  that  is  a  better  fit  to 
the  data,  particularly  for  the  trouble  spots  corresponding  to  stations  6-10  and  36- 
42,  show  that  this  model  is  more  robust,  and  improve  our  understanding  of  the 
subduction  system  beneath  Central  Washington. 

7.  Conclusion 
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We  have  generated  2-D  resistivity  models  of  the  Cascadia  subduction  zone 
structure  beneath  central  Washington  State  using  electric  and  magnetic  field  data 
collected  along  a  dense  profile  of  60  wideband  and  20  long  period 
magnetotelluric  stations.  The  resulting  images  show  a  conductive  feature  near 
the  slab  at  a  depth  of  ~40  km  that  is  consistent  with  a  fluid  release  due  to 
eclogitization  reactions  in  the  upper  crust  that  have  also  been  detected  by  a 
collocated  seismic  study.  A  second,  more  intense  conductive  feature  is  identified 
near  the  slab  at  a  depth  of  ~80  km,  consistent  with  a  fluid/melt  phase  associated 
with  the  dehydration  of  chlorite  and  serpentine  in  the  upper  mantle  and  the 
transition  of  metastable  gabbros  to  eclogite  in  the  lower  crust  of  the  descending 
slab.  This  conductor  extends  upward  to  depths  of  less  than  30  km,  where  it 
connects  with  an  intense  conductor  that  is  consistent  with  a  magma  chamber 
below  a  strong  resistive  feature  consistent  with  the  location  of  Mt.  Rainier. 

High  resistivities  are  found  beneath  the  Western  Cascades,  dropping  off 
quickly  beneath  the  more  recently  active  High  Cascades  to  the  east.  Three 
additional  shallow  conductive  features  correspond  with  sedimentary  basins 
along  the  transect. 

While  sensitivity  tests  and  comparison  with  previous  seismic  work  give  us 
confidence  in  the  robustness  and  accuracy  of  the  primary  features  in  the  model, 
we  also  have  ample  evidence  to  suggest  that  some  of  the  details  might  be  better 
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resolved.  These  issues  will  be  addressed  from  both  a  theoretical  and  practical 
standpoint  in  the  following  chapter. 
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4.  Building  a  better  MT  model:  An  investigation  into 
the  influence  of  the  assumption  of  smoothness  in 
magnetotelluric  inversion,  and  application  to  the 
Cascadia  subduction  system  in  central  Washington 


Abstract 

In  this  chapter  we  begin  by  exploring  the  assumption  of  smoothness  built  into 
to  standard  magnetotelluric  modeling,  and  the  ramifications  of  making  that 
assumption.  Through  a  series  of  models,  we  show  that  this  assumption  can 
make  it  difficult  to  resolve  sharp  gradients  in  resistivity.  We  demonstrate  that 
when  such  features  occur  in  the  true  model,  that  tell-tale  signs  in  the  inversion 
model  such  the  retreat  of  the  resistive  feature  give  us  insight  into  the  nature  of 
the  true  model,  and  that  by  using  a  combination  of  imposed  features  and  tear 
zones,  where  we  collapse  the  smoothness  assumption  locally,  we  can  find  a  more 
accurate  inverted  model.  We  also  show  that  the  assumption  of  smoothness 
inhibits  our  ability  to  resolve  the  complexity  of  a  complicated  conductive  feature, 
and  again,  we  show  that  the  imposition  of  a  tear  zone  can  help  with  this.  We 
apply  the  lessons  of  the  chapter  to  the  CAFE  data  in  order  to  generate  a  model 
that  more  clearly  shows  the  nature  of  the  conductive  features  in  the  subduction 
zone,  particularly  that  the  fluid  and  melt  phase  rises  directly  from  the  slab  at 
depths  of  80-95  km,  consistent  with  the  depths  for  dehydration  of  chlorite 
serpentine  and  the  eclogitization  of  metastable  basalts.  Assessment  of  this 
augmented  model  demonstrates  that  it  is  more  robust  than  the  model  derived  in 
Chapter  3. 
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1.  Overview 


Geophysical  models  are  never  entirely  correct.  Measurements  of 
geophysical  data  are  inevitably  inexact,  and  station  coverage  woefully 
incomplete.  Issues  of  sensitivity  and  resolution  place  boundaries  on  what  is 
detectable,  and  simplifying  assumptions  are  built  into  the  inversion 
mathematics.  While  inversion  methods  allow  us  to  create  models  that  fit  the 
data  as  closely  as  desired,  there  is  a  tradeoff  between  resolving  real  structure 
and  generating  false  structure  from  noise  that  tips  decidedly  in  favor  of  the 
latter  as  the  misfit  becomes  small.  For  any  reasonable  selection  of  misfit, 
there  are  an  infinite  number  of  solutions  that  fit  the  model  equally  well 
(Parker,  1994). 

The  selection  of  a  specific  model  is  necessarily  a  subjective  process 
involving  the  incorporation  of  a  priori  information  and/or  assumptions  about 
the  nature  of  the  distribution  of  the  measured  parameter  in  the  subsurface. 
This  subjectivity  demands  that  we  approach  each  model  critically, 
understanding  the  assumptions  that  inhabit  each  one,  and  how  they  affect  the 
model.  We  need  to  be  able  to  evaluate  the  accuracy  of  a  given  model,  and 
consider  ways  in  which  the  model  might  be  improved.  In  MT,  the  non-linear 
relationship  between  the  data  and  the  model  often  makes  it  more  difficult  to 
carry  out  these  evaluations,  but  a  fairly  standard  set  of  methods  has  arisen 
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that  MT  practitioners  utilize  to  evaluate  the  robustness  of  a  given  feature  and 
test  the  resolution  of  the  model,  and  measure  the  sensitivity  of  the  model  to 
perturbations. 

Non-linearity  also  implies  that  even  were  the  data  error-free,  it  would  not 
be  possible  for  the  model  parameters  to  be  equally  well  determined,  requiring 
careful  study  of  the  sensitivity  of  the  data  (e.g.,  Schwalenberg  et  al.,  2002). 
Pseudo-sections,  which  show  a  profile  of  apparent  resistivity  and  phase 
interpolated  over  both  period  and  distance,  are  a  useful  tool  for  exploring 
data  sensitivity.  Software  packages  such  as  WinGLink  allow  for  model 
response  derived  synthetic  pseudo-sections,  which  can  then  be  differenced 
from  pseudo-section  generated  from  the  data.  The  resulting  residual  pseudo¬ 
sections  highlight  regions  for  which  the  forward  model  does  not  closely 
reproduce  the  data.  Some  examples  of  this  are  included  in  Appendix  2a. 

It  has  also  been  possible  to  use  the  elements  of  the  sensitivity  matrix 
directly  as  a  measure  of  how  small  model  distortions  affect  the  data 
(Schwalenberg  et  al.,  2002,  Baba  et  al.,  2006),  but  as  this  is  a  strictly  linear 
approach  it  only  has  value  in  very  close  proximity  to  the  derived  model. 
Furthermore,  full  calculation  of  this  sparse  matrix  at  every  step  is  very 
computer-intensive  and  for  modern  inversion  methods  such  as  non-linear 
conjugate  gradients  the  sensitivity  matrix  is  only  considered  within  matrix- 
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vector  products,  and  neither  fully  calculated  nor  retained  at  every  step  (Rodi 
and  Mackie,  2001). 

One  of  the  most  important  tools  used  by  MT  investigators  to  explore  a 
given  model  is  manual  editing  of  the  model,  which  can  be  implemented  in 
either  a  forward  or  inverse  sense.  As  a  forward  modeling  tool,  manual 
editing  serves  to  explore  the  range  of  models  that  are  allowed  by  the  data  by 
varying  different  parameters  such  as  resistivity  or  extent  of  a  feature  in  a 
systematic  way  and  observing  the  effect  on  the  misfit  between  the  data  and 
model.  For  example,  in  a  subduction  zone  setting  with  a  robust  volcanic 
front,  we  might  suspect  a  conductive  fluid  release  from  the  slab  beneath  the 
volcanic  arc,  even  if  it  is  not  apparent  from  the  initial  inversion  model.  We 
can  decrease  the  resistivity  in  the  region  above  the  slab  until  the  misfit 
becomes  worse,  thereby  establishing  a  limit  for  the  lowest  resistivity  allowed 
by  the  data.  This  has  been  done  in  Oregon  to  show  that  a  highly  conductive 
feature  is  allowed  (but  not  required)  by  the  data  (Evans  et  al.,  submitted),  and 
in  Costa  Rica,  reducing  the  resistivity  was  shown  to  even  slightly  improve  the 
fit  between  the  data  and  model  (Worzewski  et  al.,  2010). 

This  method  can  also  be  used  to  find  the  allowable  extent  of  a  given 
feature.  In  southern  Tibet,  the  shallowest  allowable  bottom  of  a  poorly 
constrained  conductor  was  located  by  imposing  a  resistive  halfspace  at  depth. 
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and  raising  the  depth  of  the  halfspace  until  the  misfit  worsened,  enabling  the 
investigators  to  constrain  a  maximum  value  for  resistivity  as  well  (Li  et  al., 
2003).  Beneath  the  East  Pacific  Rise,  the  depth  to  a  conductive  layer  and  the 
lateral  extent  of  a  conductive  plume  feature  were  constrained  conjunctively 
by  evaluating  the  misfit  while  varying  these  parameters  (Baba  et  al.,  2006). 
While  widely  used  and  undeniably  useful,  these  forward  modeling  tools  are 
limited  by  the  non-linearity  of  MT  in  that  they  do  not  account  for  the 
inductive  coupling  between  features. 

As  an  inversion  tool,  conductive  and  resistive  features  can  be  removed 
from  a  model,  and  the  inversion  re-run,  either  without  constraints  or  seeking 
a  model  near  the  starting  model.  If  the  feature  recurs  during  the  inversion  or 
the  fit  is  not  reacquired,  then  the  feature  is  taken  to  be  required  by  the  data.  If 
the  feature  does  not  recur  or  comes  back  in  an  altered  sense,  and  the  fit  is 
reacquired,  we  obtain  information  about  the  range  of  models  that  the  data 
will  accept.  Extensive  application  of  this  method  to  key  features  of  an  MT 
model  is  standard  practice;  one  rather  detailed  example  can  be  found  in  a 
study  of  the  structure  of  the  lithosphere  in  NE  Botswana  (Miensopust  et  al., 
2011). 

This  method  can  also  be  applied  to  incorporate  a  priori  information  into  an 
inversion  sequence.  A  feature  can  be  built  into  the  starting  model  and  fixed. 


110 


as  is  standard  practice  when  an  ocean  is  present  in  a  model,  or  built  into  the 
starting  model  and  allowed  to  evolve  with  the  inversion.  In  the  latter  case, 
the  idea  is  to  constrain  the  starting  "neighborhood"  to  give  the  inversion  a 
better  opportunity  to  operate  nearer  the  global  minimum.  For  example,  it  is 
well  known  that  MT  can  be  rather  insensitive  to  the  resistivity  of  a  highly 
resistive  and  extensive  feature  such  as  a  subducting  slab,  so  by  including  a 
resistive  feature  in  the  presumed  location  of  the  slab,  it  is  possible  to  constrain 
the  inversion  sequence  to  search  for  solutions  that  include  that  feature.  This 
can  illuminate  other  model  features  that  might  otherwise  have  remained 
undetected  or  poorly  constrained  (Brasse  et  al.,  2009,  Matsuno  et  al.,  2010, 
Evans  et  al.,  submitted). 

The  combination  of  manual  editing  and  re-inversion  provides  a  powerful 
and  non-linear  technique  for  exploring  the  range  of  acceptable  models,  but  it 
does  nothing  to  challenge  the  assumptions  inherent  in  the  inversion  itself.  In 
this  chapter  we  will  examine  potential  consequences  of  the  assumption  of 
smoothness  that  is  built  into  the  inversion  process,  and  explore  some 
techniques  that  can  be  utilized  to  recognize  and  mitigate  some  of  those 
consequences.  Finally,  we  will  apply  those  lessons  to  the  data  from  Chapter 
3,  with  the  objective  of  producing  a  demonstrably  more  accurate  MT 
inversion  model. 
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All  of  the  models  in  this  chapter  were  generated  using  the  WinGLink  MT 
software  package  (Rodi  and  Mackie,  2001).  To  standardize  the  models  and 
maximize  the  applicability  to  the  CAFE  MT  data  set,  the  same  parameters 
were  used  for  all  of  the  inversions  unless  otherwise  noted.  The  profiles  are  all 
300  km  long  with  60  stations  arranged  at  roughly  5  km  distance,  and  the 
inversion  grids  are  all  76  rows  by  196  columns,  with  row  width  increasing 
with  depth.  All  inversions  solve  for  the  smoothest  model  using  the  uniform 
grid  Laplacian  operator.  Error  floors  are  set  to  five  percent  for  the  TM  mode 
and  ten  percent  for  the  TE  mode,  and  the  tipper  data  is  included.  The 
parameters  a  and  [3,  which  define  the  relationship  of  the  smoothness 
parameter  in  the  vertical  as  compared  to  the  horizontal,  and  the  way  in  which 
the  smoothness  changes  with  depth  respectively,  are  set  to  1.5  and  1.7 
respectively,  and  the  regularization  parameter  x  is  equal  to  3.3. 

2.  Theoretical  underpinnings 

The  general  2-D  MT  inverse  problem  can  be  thought  of  in  terms  of 
attempting  to  solve  for  m  in  the  equation 

d  =  F(m)  +  e  (1) 
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where  the  components  of  the  N-dimensional  column  vector  d  represent  the 


observed  magnetic  and  electric  fields  (in  the  form  of  amplitude  and  phase  of 
apparent  resistivity  for  a  given  station,  period,  and  polarization)  and  the 
components  of  the  M-dimensional  column  vector  m  represent  the  resistivity 
of  the  individual  elements  in  model  space.  F  is  the  transformation  from 
model  space  to  data  space  required  by  Maxwell's  equations,  and  the 
components  of  the  N-dimensional  column  vector  e  represent  the  data  errors. 

Unfortunately,  solving  equation  (1)  for  m  turns  out  to  be  impossible  for  a 
number  of  reasons,  the  first  of  which  is  that  the  chances  of  a  solution  actually 
existing  is  almost  certain  to  be  zero,  especially  if  the  column  vector  e  *  0.  The 
standard  approach  would  be  to  move  away  from  the  idea  of  an  exact 
solution,  and  instead  to  redefine  solution  to  mean  an  m  for  which  F(m) 
matches  d  as  closely  as  possible  in  some  pre-defined  way.  This  is  done  by 
posing  the  problem  as  an  optimization  problem,  where  an  objective  function 
is  minimized  subject  to  various  constraints. 

vF(m)  =  (d  -  F(m))T  V  (d  -  F(m))  (2) 
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In  a  least  squares  approach,  such  as  that  shown  in  equation  (2),  the  objective 
function  VF  is  defined  as  the  squared  L2norm  of  the  residual  vector,  subject  to 
a  positive-definite  weighting  matrix  V,  which  is  typically  related  to  the 
variance  of  the  error  vector. 

This  brings  us  to  our  second  difficulty,  namely  that  in  order  to  ensure  the 
existence  of  a  least  squares  solution  to  a  real  geological  problem,  we  must 
either  make  the  size  of  m  very  large  by  using  a  dense  sampling  grid,  or  define 
the  solution  set  so  broadly  as  to  not  be  useful.  In  either  case  we  ensure  that 
any  solution  that  exists  is  not  unique,  and  so  the  least  squares  problem 
becomes  ill-posed  and  unstable.  For  inverse  problems  involving  large  model 
spaces,  the  standard  approach  is  to  introduce  a  stabilizing  term  that  favors 
some  models  over  others  by  penalizing  some  pre-defined  property  that  is 
deemed  to  be  undesirable.  While  several  stabilizing  functions  have  been 
successfully  used  in  various  geophysical  applications,  many  magnetotelluric 
algorithms  employ  the  minimum  norm  of  the  Laplacian,  which  will  seek  the 
smoothest  model  possible  that  fits  the  data.  The  argument  is  that  by  seeking 
the  smoothest  model  (for  a  given  misfit),  only  structure  explicitly  required  by 
the  data  will  be  resolved  (de  Groot-Hedlin  and  Constable,  1990).  In  practice, 
this  stabilizing  function  rules  out  many  geologically  incomprehensible 
mathematical  solutions  and  forces  the  algorithm  to  seek  solutions  that  search 
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for  structures  with  coherent  regional  features.  The  regularized  solution  that 
is  used  to  generate  all  of  the  models  in  this  study  is  given  in  equation  (3) 
below 


W(m)  =  (d  -  F(m))T  V  (d  -  F(m))  +  AmTLTLm  (3) 

where  the  matrix  L  is  a  simple,  second-difference  operator  such  that  Lm 
approximates  the  log  of  resistivity  for  a  uniform  grid  of  model  blocks  (Rodi 
and  Mackie,  2001).  The  regularization  parameter  A  is  a  positive  number  that 
determines  globally  the  weight  given  to  the  stabilizing  functional.  The 
selection  of  A  is  chosen  so  as  to  balance  model  roughness  with  data  misfit.  It 
is  generally  chosen  either  with  a  specific  target  misfit  or  at  such  a  place  on  the 
misfit  vs.  roughness  curve  that  increasing  A  results  in  a  small  increase  in 
model  roughness  relative  to  the  decrease  in  misfit,  and  decreasing  it  results  in 
a  small  decrease  in  misfit  relative  to  the  increase  in  roughness. 

Finally,  the  non-linearity  of  the  problem  requires  that  the  solution  (the 
minimization  of  the  objective  function  VF)  be  approached  iteratively.  The 
non-linear  conjugate  gradient  (NLCG)  method  used  in  this  study  is  a 
modification  of  the  steepest  descent  method,  where  the  search  direction  in  an 
N-dimensional  topography  is  defined  as  the  direction  along  which  the  value 
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of  the  objective  function  decreases  most  rapidly  and  proceeds  until  it  is 
minimized  along  that  line.  At  that  point,  it  searches  for  the  direction  of 
maximum  decrease  again.  For  NLCG,  the  search  direction  is  a  linear 
combination  of  the  line  of  steepest  descent  and  the  previous  search  direction. 
Going  "downhill"  along  these  planes  not  only  helps  to  avoid  small  local 
minima  (Zhdanov,  1993),  but  greatly  reduces  the  calculation  effort  by 
retaining  directional  information  from  previous  iterations  within  each  new 
search  direction. 

3.  A  closer  look  at  the  stabilizing  functional 

The  effect  of  the  stabilizing  functional  shown  in  equation  (3)  is  to  penalize 
models  in  which  the  conductivity  changes  sharply.  Not  only  does  this 
stabilize  the  inversion  and  allow  us  to  find  a  "best"  model  for  a  given  misfit, 
but  it  prevents  an  entire  class  of  geologically  unsound  models  in  which 
anomalous  conductivities  that  are  purely  artifacts  of  the  inversion  process 
populate  model  space.  While  this  is  clearly  desirable,  it  is  also  true  that  sharp 
boundaries  do  exist  in  nature,  and  smooth  algorithms  may  not  sufficiently 
resolve  these  boundaries. 
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In  order  to  address  these  difficulties,  several  algorithms  have  been 
developed  that  allow  for  sharp  boundaries  by  making  different  global 
assumptions  about  the  distribution  of  the  modeled  parameter  in  the 
subsurface,  for  example  by  favoring  models  made  up  of  layers  (Smith  et  al., 
1999,  de  Groot-Hedlin  and  Constable,  2004)  or  containing  primarily  blocky 
structure  (Mehanee  and  Zhdanov,  2002).  Of  course,  each  of  these  global 
algorithms  present  trade-offs,  and  are  effective  largely  to  the  extent  that  the 
underlying  assumptions  are  accurate  throughout  the  modeled  space.  While  a 
priori  knowledge  regarding  the  location  and  nature  of  sharp  discontinuities 
can  be  an  extremely  useful  implement  in  resolving  the  associated  structure, 
its  utility  in  resolving  features  in  a  sufficiently  complicated  subsurface  is 
dependent  on  the  extent  to  which  it  can  be  effectively  applied  locally. 

Inversion  packages  such  as  WinGLink  (Rodi  and  Mackie,  2001)  and 
Occam  (de  Groot-Hedlin  and  Constable  1990)  allow  for  the  elimination  of  the 
penalty  for  roughness  along  presumed  discontinuities,  and  while  this  tool  has 
recently  been  used  to  good  effect  (Matsuno  et  al.,  2010,  Miensopust  et  al., 

2011,  Evans  et  al.,  submitted),  we  are  cautioned  that  incorrect  placement  of 
such  boundaries  may  produce  misleading  results  (de  Groot-Hedlin  and 
Constable,  1990).  While  this  is  certainly  true,  we  must  also  recognize  that  the 
use  of  such  tools  is  simply  a  matter  of  trading  one  set  of  assumptions  for 
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another  and  we  should  be  cautious  in  any  event.  We  should  also  be  asking  to 
what  extent  the  assumption  of  smoothness  can  return  misleading  results. 

Figure  1  shows  the  results  of  an  experiment  designed  to  provide  one 
answer  to  that  question.  The  forward  model  (fig  la)  is  a  100  Cm  halfspace  in 
which  a  10,000  Qm  resistive  block  and  a  1  Cm  conductive  block  are  placed  in 
contact  with  one  another  from  depths  of  60-115  km.  Each  block  is  60  km 
wide.  The  forward  response  to  this  conductivity  profile  was  generated  and 
saved  as  synthetic  data  to  an  array  of  60  stations  located  at  roughly  5  km 
intervals  at  the  surface. 

We  then  generated  inversions  from  the  synthetic  data  using  the  WinGLink 
software  package.  Instead  of  using  a  halfspace  starting  model  however,  we 
used  a  perfect  copy  of  the  forward  model,  such  that  initially,  the  (F(m)-d) 
term  was  equal  to  zero.  Without  the  stabilizing  functional,  the  objective 
function  would  also  have  been  equal  to  zero,  and  the  inversion  would  have 
stopped  immediately  with  the  correct  solution.  With  the  stabilizing 
functional,  despite  the  fact  that  the  misfit  for  the  first  iteration  was  zero,  the 
objective  function  was  not  minimized  and  the  inversion  proceeded.  The 
implication  is  that  the  differences  between  the  forward/starting  model  (Fig 
la)  and  the  inversion  model  (Fig.  lb)  can  be  entirely  attributed  to  the  effect  of 
the  smoothing  stabilizing  functional  on  the  inversion. 
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The  conductor  has  for  the  most  part  been  retained  in  the  correct  location, 
with  an  accurate  resistivity  of  lQm  near  the  center.  The  corner  detail  of  the 
original  feature  have  been  lost,  as  has  its  uniformity,  as  the  intensity  of  the 
conductive  feature  in  the  inversion  grades  out  in  roughly  concentric  circles. 


a. 
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Figure  1:  Inversion  model  demonstrating  the  effect  of  a  smoothing  stabilizing 
functional  for  a  simple  two  block  model.  Figure  la  shows  the  forward  model, 
which  was  also  used  as  the  starting  model,  and  Figure  lb  shows  the  inversion. 
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The  resistive  block  has  fared  less  well,  with  the  highest  resistivity  more  than  a 
full  order  of  magnitude  too  low  at  600  Qm  and  displaced  nearly  15  km  from  the 
original  center,  fully  55  km  from  the  most  conductive  value  of  1  Qm  despite 
sharing  a  common  boundary  in  the  forward  and  starting  models.  The  resistor 
appears  stretched  out  in  the  vertical,  wrapping  around  the  conductive  structure 
at  a  nearly  constant  distance,  and  a  phantom  resistor  of  150-200  Qm  has 
appeared  on  the  opposite  side  of  the  conductor. 

None  of  this  is  intended  to  suggest  that  the  incorporation  of  the  smoothing 
stabilizing  functional  is  misguided  in  any  way.  In  the  absence  of  a  priori 
information,  it  has  been  the  most  effective  way  to  stabilize  the  objective  function 
and  ensure  that  features  in  the  inversion  are  required  by  the  data.  The  starting 
model  in  our  example  was  specifically  selected  to  highlight  the  limitations  of  a 
smooth  model  assumption,  and  a  general  interpretation  of  the  final  image  would 
still  provide  good  insights  into  the  conductivity  structure  of  the  subsurface. 
Instead,  we  strive  to  develop  a  better  understanding  of  the  tradeoffs  inherent  in 
this  assumption,  and  to  explore  methods  that  we  can  employ  to  improve  the  final 
model,  particularly  as  they  apply  to  a  subduction  setting. 

4.  Exploration 
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In  this  section,  we  explore  the  effects  of  incorporating  a  priori  information  into 
the  inversion  of  synthetic  data  generated  from  known  models.  We  looked  at 
more  than  200  full  inversions  from  dozens  of  models  ranging  in  complexity  from 
simple  layers  and  blocks  to  simulated  subduction  zone  models.  Two  techniques 
proved  to  be  particularly  useful  both  to  improve  the  accuracy  of  the  inversion 
models,  and  to  provide  additional  information  about  the  subsurface  to  assist 
interpretation. 

The  first  of  these  is  to  impose  high  resistivity  in  the  initial  models  at  locations 
where  it  would  be  expected,  such  as  the  location  of  a  subducting  slab,  but 
leaving  the  feature  free  to  evolve  with  the  model.  Conceptually,  the  idea  is  to 
start  the  inversion  in  a  topographical  neighborhood  that  is  closer  to  the  global 
minimum  for  the  objective  function  than  the  neighborhood  representing  a  simple 
halfspace.  In  practice,  the  presence  of  the  resistor  tends  to  help  illuminate  nearby 
conductors.  Additionally,  as  smoothness  is  imposed  on  the  inversion,  the 
boundary  of  an  imposed  resistor  will  appear  to  retreat  from  a  nearby  conductor 
in  order  to  minimize  model  roughness;  the  extent  of  this  retreat  is  positively 
correlated  to  the  intensity  and  proximity  of  the  conductor.  Conversely,  the 
imposition  of  a  conductive  feature  in  starting  models  did  not  seem  to  improve 
the  inversion. 
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The  second  technique  was  the  imposition  of  a  zone  in  which  the  smoothness 
constraint  was  collapsed  locally,  hereafter  referred  to  as  a  tear  zone.  These  tear 
zones  were  primarily  imposed  in  places  where  we  suspected  a  rapid  change  in 
conductivity,  such  as  at  the  top  of  a  subducting  slab  where  dehydration  reactions 
are  occurring  and  fluids  are  being  released.  Tear  zones  also  proved  useful  in 
illuminating  the  structure  of  a  conductor,  as  the  imposed  smoothness  tends  to 
obfuscate  detail  and  ensure  that  the  most  conductive  region  appears  near  the 
center. 

We  did  not  succeed  in  developing  any  method  that  works  in  every  situation; 
presumably  the  non-linear  nature  of  MT  precludes  this  possibility.  What  we 
have  instead  are  a  set  of  tools  that  can  be  used  first  to  develop  an  intuition  about 
a  given  problem,  and  then  to  work  towards  an  improved  inversion  model. 
Finally,  we  can  evaluate  the  models  and  demonstrate  which  of  them  likely  most 
accurately  reflects  the  conductivity  structure  of  the  subsurface. 

The  first  set  of  inversion  models  is  shown  in  figure  2.  The  forward  model  (fig. 
2a)  includes  a  9000  Qm  dipping  resistive  layer  intended  to  represent  a  slab,  with 
two  relatively  thin  3  Qm  conductors  at  the  top  of  the  slab  that  correspond  with 
where  we  might  expect  fluids  released  during  dehydration  reactions  to  occur.  A 
0.33  Qm  conductive  feature  representing  the  ocean  is  also  present  in  the  upper 
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left  corner  of  the  model.  This  conductor  is  locked  in  for  all  of  the  inverted 
models  just  as  would  be  for  an  inversion  of  real  data  in  a  coastal  MT  study. 


Figure  2:  Inversion  model  for  a  subducting  slab  with  thin  conductors  at  depths 
corresponding  to  fluid  release  from  dehydration  reactions  for  Cascadia.  Fig.  2a 
shows  the  forward  model.  Fig.  2b  shows  the  inversion  model  from  a  halfspace. 
Fig.  2c  shows  the  inversion  model  with  the  resistive  slab  imposed  into  the 
starting  model.  Fig.  2d  shows  the  inversion  model  with  the  imposed  resistive 
slab  and  a  tear  zone  incorporated  into  the  slab. 
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The  halfspace  inversion  model  (fig.  2b)  suffers  from  some  of  the  same 
limitations  that  we  observed  in  the  two  block  model.  Both  conductors  are 
resolved  in  essentially  their  correct  locations,  but  much  of  the  detail  of  their 
shapes  has  not  been  captured.  The  dipping  resistive  layer  is  essentially  non¬ 
existent  downstream  of  the  shallow  conductor.  This  smoothing  away  of  the 
resistive  feature  also  results  in  overestimating  the  resistivity  of  the  two 
conductors;  the  lowest  resistivity  in  the  shallow  conductor  is  10  Qm  and  in  the 
deeper  conductor  24  Qm. 

Figure  2c  shows  the  model  for  the  inversion  in  which  the  starting  model 
included  the  dipping  resistor.  While  the  inversion  misfit  tends  not  to  be  terribly 
sensitive  to  changes  to  the  resistive  feature,  the  inclusion  of  the  resistor  in  the 
starting  model  helps  to  illuminate  the  details  of  the  overlying  conductors.  The 
location  of  both  conductors  can  now  clearly  be  seen  to  be  at  the  top  of  the 
dipping  resistor.  Additionally,  both  are  better  resolved  in  space  and  suffer  from 
less  smearing.  Their  resistivities  are  also  approaching  their  correct  values,  with 
the  resistivity  of  the  shallow  conductor  at  7  Qm  and  that  of  the  deep  conductor 
17  Qm.  The  dipping  resistive  feature  is  also  still  present  in  the  model,  but  with  a 
highest  resistivity  of  only  1000  Qm.  Additionally,  it  has  retreated  away  from  the 
conductors  just  as  we  saw  for  the  two  block  model;  we  have  come  to  recognize 
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this  as  a  reliable  indicator  that  the  distance  between  the  conductors  and  resistor 


in  the  starting  model  (or  by  extension  the  real  subsurface)  is  not  very  large. 

Our  final  inversion  model  for  this  series  is  shown  in  figure  2d.  For  this 
inversion,  we  imposed  the  dipping  resistor  in  place,  and  also  placed  a  tear  zone 
covering  the  resistor,  encouraging  the  inversion  to  seek  solutions  where  the 
conductivity  can  change  abruptly  at  its  edges.  This  model  has  a  very  high  level 
of  accuracy.  The  shallow  conductor  is  located  in  the  correct  place  at  the  top  of 
the  slab,  with  a  resistivity  of  as  little  as  3  Cm.  The  deeper  conductor  is  also 
correctly  located,  with  resistivities  down  to  10  Qm. 

The  imposition  of  the  tear  zone  at  the  upper  boundary  of  the  dipping  resistor 
changes  the  assumption  implicit  in  the  inversion  from  seeking  models  with 
smooth  transitions  to  those  with  sharp  boundaries.  From  forward  models  for 
which  the  conductors  share  a  boundary  with  the  dipping  resistor,  this  leads 
towards  improved  accuracy  in  the  inversion,  but  we  should  also  explore  models 
in  which  the  boundary  is  not  shared  to  understand  what  potential  the  imposition 
of  a  tear  zone  has  to  lead  us  astray. 

The  series  of  inversions  in  figure  3  is  designed  to  help  answer  this  question. 
The  forward  model  (fig.  3a)  consists  of  a  9000  Qm  dipping  resistive  feature  and 
0.33  Qm  ocean  identical  to  those  in  the  previous  example.  There  is  also  a  pair  of 
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3  Qm  conductors,  but  in  this  case,  they  have  been  separated  from  the  dipping 


resistor. 


Figure  3:  Inversion  model  for  a  subducting  slab  with  conductors  separated  from 
the  top  of  the  slab.  Fig.  3a  shows  the  forward  model.  Fig.  3b  shows  the  inversion 
model  from  a  halfspace.  Fig.  3c  shows  the  inversion  model  with  the  resistive  slab 
imposed  into  the  starting  model.  Fig.  3d  shows  the  inversion  model  with  the 
imposed  resistive  slab  and  a  tear  zone  incorporated  into  the  slab. 
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The  halfspace  model  (fig.  3b)  does  reasonably  well  recovering  all  of  the 
features  in  the  inversion.  Both  conductors  are  present  and  in  the  correct  location, 
and  while  the  deeper  conductor  has  faded  to  7  Qm,  the  shallow  one  has 
maintained  its  original  resistivity  of  3  Qm.  Even  the  dipping  resistor  is  resolved 
here  to  some  degree,  with  resistivities  of  several  hundred  Qm  extending  beyond 
depths  of  80  km. 

Figure  3c  shows  an  inversion  model  in  which  the  dipping  resistor  was 
imposed  into  the  starting  model.  Both  conductors  are  very  well  resolved,  with 
the  deeper  of  the  two  reaching  a  resistivity  of  4  Qm.  The  dipping  resistive 
feature  remains  largely  intact  as  well,  with  resistivities  of  more  than  1000  Qm 
extending  to  depths  of  more  than  120  km.  The  resistor  also  shows  very  little  sign 
of  retreat  due  to  smoothing  out  in  the  vicinity  of  a  strong  conductor.  This  is  a 
good  indication  that  the  conductors  do  not  share  a  boundary  with  the  dipping 
resistive  feature,  and  counter  indicative  to  installing  a  tear  zone  at  the  top  of  the 
conductor. 

Nevertheless,  figure  3d  shows  an  inversion  model  for  which  a  tear  zone  was 
imposed  over  the  dipping  resistor.  The  shallow  conductor  maintains  a  low 
resistivity  value  of  3  Qm  in  the  correct  location,  but  also  appears  to  be  connected 
to  the  slab  feature,  which  we  know  not  to  be  the  case.  The  situation  for  the 
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deeper  conductor  is  even  more  problematic,  with  resistivities  as  low  as  6  Cm  that 
appear  to  be  in  direct  contact  with  the  dipping  resistor. 

The  two  examples  shown  are  fairly  clear  cut,  but  it  is  obviously  possible  for 
intermediate  cases  to  exist  for  which  determining  whether  or  not  to  use  a  tear 
zone  would  be  more  challenging.  While  the  careful  imposition  of  resistors  in  the 
starting  model  can  simplify  the  task  by  constraining  the  location  of  the 
conductive  features,  there  are  other  tools  that  we  can  use  to  evaluate  the  final 
inverted  models.  We  will  explore  this  more  fully  in  section  5. 

For  the  next  set  of  inversion  models  (Fig.  4),  we  move  on  to  another 
application  of  the  tear  zone.  Figure  4a  shows  a  forward  model  consisting  of  a 
9000  Cm  resistive  layer,  and  a  conductive  feature  made  up  of  two  1  Cm  blocks 
that  overlap  at  one  corner. 

The  halfspace  inversion  model  is  shown  in  fig.  4b.  The  detail  in  the  original 
conductor  has  been  lost  to  the  smoothness  constraints,  replaced  by  a  single 
conductor  with  a  core  of  resistivity  1  Cm  in  the  upper  right  quadrant, 
surrounded  by  concentric  elliptical  shapes  of  increasing  resistivity,  with  the 
whole  occupying  roughly  the  same  area  as  the  original  conductor. 

The  inversion  model  in  fig.  4c  uses  a  starting  model  with  a  tear  zone  imposed 
over  the  area  occupied  by  the  conductor  in  the  halfspace  model,  resulting  in  a 
partial  recovery  of  the  detail  of  the  original  feature.  The  conductor  in  this  model 
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clearly  exhibits  the  dual  centers  of  the  forward  model,  one  in  the  upper  right  and 


the  other  in  the  lower  left,  both  with  resistivities  in  the  1-2  Qm  range. 
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Figure  4:  Inversion  models  demonstrating  the  loss  of  detail  for  a  conductive 
feature  due  to  the  smoothness  constraint,  and  the  use  of  tear  zones  to  recover 
some  of  the  detail.  Fig.  4a  shows  the  forward  model.  Fig.  4b  shows  a  halfspace 
inversion  model.  Fig.  4c  shows  an  inversion  model  in  which  a  tear  zone  has  been 
imposed  over  the  boundaries  of  the  conductor. 
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While  the  smoothness  constraints  tend  to  affect  most  conductive  features,  it  is 
not  always  possible  to  recover  detail  using  tear  zones.  One  set  of  inversions  for 
which  this  strategy  is  ineffective  is  shown  in  figure  5.  The  forward  model  in  fig. 
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5a  shows  a  conductive  column  with  conductivities  that  increase  with  depth.  The 
bottom  layer  has  a  resistivity  of  1  Qm,  with  layers  of  3  Qm,  8  Qm,  and  14  Qm 


Figure  5:  Inversion  models  demonstrating  the  loss  of  detail  in  a  conductor  with  a 
gradient.  Fig.  5a  shows  the  forward  model.  Fig.  5b  shows  the  halfspace 
inversion.  Fig.  5c  shows  an  inversion  with  a  tear  zone  imposed  over  the  region 
of  the  conductor,  and  Fig.  5d  an  inversion  model  with  both  an  imposed  resistive 
layer  and  tear  zone  over  the  resistive  layer. 
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respectively.  The  base  of  the  bottom  layer  is  in  contact  with  a  thick  resistive 
layer  with  a  resistivity  of  9000  Cm. 

We  can  see  in  the  halfspace  inversion  (fig.  5b)  that  the  conductive  gradient  in 
the  forward  model  has  been  replaced  by  the  characteristic  conductive  feature 
that  we  have  now  seen  many  times,  a  concentric  series  of  conductivities  that 
decrease  outward  from  the  center,  this  time  with  resistivities  as  low  as  4  Cm. 

The  underlying  resistive  feature  is  weakly  in  evidence,  with  some  suggestion  of 
the  retreating  character  that  we  have  come  to  identify  with  a  conductor  in  close 
proximity.  This  can  be  easily  confirmed  by  running  an  inversion  with  the 
imposition  of  a  resistive  layer  (which  was  conducted  but  the  results  are  not 
shown  here  as  they  demonstrate  nothing  new). 

Our  attempt  to  recapture  the  detail  of  the  original  conductor  by  imposing  a 
tear  zone  on  the  area  of  the  conductor  is  shown  in  fig.  5c.  In  this  case,  we  are 
unable  to  improve  the  image,  as  the  influence  of  the  nearby  resistor  is  too  strong. 
If  the  conductor  were  isolated  in  the  halfspace,  we  would  be  able  to  improve  the 
resolution  to  some  degree,  but  with  the  resistive  layer  directly  below,  imposing 
the  tear  on  the  conductor  simply  allows  it  to  contract  in  the  upper  region  of  the 
tear  zone. 

Since  the  imposed  resistive  layer  retreats  during  the  inversion,  we  expect  that 
the  conductive  feature  is  very  close,  leading  us  to  try  imposing  a  tear  zone  over 
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the  resistor.  The  results  of  this  inversion  are  shown  in  fig.  5d.  While  this  is  the 
most  accurate  of  our  inversion  models  for  this  set,  the  details  of  the  conductor 
are  essentially  the  same  as  they  would  have  been  for  a  conductor  of  uniform 
resistivity  in  the  forward  model,  so  we  cannot  claim  that  any  detail  has  been 
salvaged. 

The  preliminary  models  for  the  CAFE  data  set  in  Chapter  3  show  a  moderate 
conductor  near  the  slab  at  a  depth  of  ~40  km,  and  a  strong  conductor  which 
appears  to  rise  from  near  the  slab  at  a  depth  of  ~80  km.  This  is  in  contrast  with 
the  results  of  another  study  in  Oregon  along  the  EMSLAB  line  (Evans  et  al., 
submitted)  which  shows  a  strong  conductor  rising  from  the  slab  at  the  shallower 
depth,  and  extending  horizontally  at  ~25  km  depth  across  the  profile  at  least  as 
far  as  the  volcanic  arc.  The  deeper  conductor  near  80  km  depth  does  not  show 
up  in  the  standard  inversions  of  that  model,  but  if  such  a  conductor  is  imposed 
into  the  starting  model,  it  is  retained  in  the  inversion  and  an  acceptable  fit  to  the 
data  is  achieved.  Inversions  of  the  tipper  data  (comparing  only  the  vertical  and 
horizontal  magnetic  fields)  alone  also  suggest  that  such  a  conductor  does  exist. 
The  question  becomes  whether  the  deeper  conductor  is  much  weaker  in  the 
Oregon  line,  or  whether  the  presence  of  the  extended  shallower  conductor 
renders  it  less  detectable. 
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The  inversion  set  in  figure  6,  while  in  keeping  with  our  general  theme,  is  part 
of  a  series  of  models  that  were  generated  specifically  to  address  this  difference 


Figure  6:  Inversion  models  demonstrating  the  effects  of  a  strong  shallow 
conductor  that  extends  horizontally  in  the  direction  of  the  volcanic  arc.  Fig.  6a 
shows  the  forward  model.  Fig.  6b  shows  the  inversion  model  from  a  halfspace. 
Fig.  6c  shows  the  inversion  model  with  a  dipping  resistor  imposed  on  the  starting 
model,  and  Fig.  6d  includes  both  the  dipping  resistor  and  an  imposed  tear  zone. 
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between  the  CAFE  and  EMSLAB  lines.  Figure  6a  shows  the  forward  model,  with 
conductor  is  actually  quite  well  resolved,  and  the  telltale  retreat  of  the  resistive 
layer  is  in  evidence  as  well.  The  deep  conductor  has  resolved  to  the  now  familiar 
concentric  shape,  with  a  lowest  resistivity  of  14  Qm. 

Imposing  the  dipping  resistor  on  the  starting  model  does  not  improve  the 
results  very  much,  although  it  does  provide  some  useful  information.  The 
shallow  conductor  loses  some  of  its  shape  as  the  imposed  smoothness  causes  the 
inversion  to  try  to  put  distance  between  it  and  the  resistor.  The  deeper 
conductor  has  a  lowest  resistivity  now  of  11  Qm,  but  its  impact  on  the  imposed 
resistor  should  cause  us  to  consider  a  more  conductive  value.  Additionally,  the 
retreat  of  the  resistor  from  both  conductors  suggests  that  we  might  do  well  to 
impose  a  tear  zone  at  the  top  of  the  dipping  resistive  feature. 

The  starting  model  for  the  inversion  image  in  figure  6d  included  both  a 
dipping  resistor  and  a  tear  along  the  top  layer  of  the  resistor.  The  model  is  quite 
accurate  with  respect  to  the  forward  model.  The  shallow  conductor  is  very  well 
resolved,  and  the  deeper  conductor  has  resistivities  as  low  as  6  Qm,  although  the 
details  of  the  shape  of  the  conductor  have  not  been  captured.  It  is  clear  from  this 
inversion  series  that  if  there  is  a  deeper  conductor  for  the  EMSLAB  transect,  that 
it  is  not  nearly  as  pronounced  as  the  deep  conductor  in  the  CAFE  line. 
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5.  Application  to  the  CAFE  data  set 


The  imposition  of  resistive  features  and  tear  zones  into  a  starting  model  have  been 
shown  to  be  effective  techniques  to  improve  the  quality  and  resolution  of  an  inversion, 
but  optimally  require  a  priori  data  concerning  the  location  of  the  top  of  the  resistive 
feature.  For  the  CAFE  study,  we  analyzed  data  collected  at  dense  arrays  of  roughly 
collocated  passive  seismic  and  magnetotelluric  stations.  Both  sets  of  data  were 
analyzed  independently,  employing  Generalized  Radon  Transform  (GRT)  migration  for 
the  seismic  data  (see  Chapter  2),  and  Non  Linear  Conjugate  Gradient  (NLCG)  inversion 
for  the  MT  (see  Chapter  3). 

There  are  a  number  of  advantages  in  investigating  collocated  seismic  and  MT  data. 
First  of  all,  because  the  two  techniques  are  sensitive  to  entirely  different  sets  of  rock 
properties  (Hellfrich,  2003),  each  is  likely  to  be  able  to  provide  information  pertaining  to 
the  structure  of  the  subsurface  that  the  other  cannot.  Secondly,  because  velocity  and 
conductivity  are  often  in  fact  correlated  (Stanley  et  al.,  1990,  Jones  et  al.,  2009)  either 
directly  or  simply  in  the  fact  that  at  distinct  geological  boundaries  both  parameters  are 
likely  to  change  (Gallardo,  2007,  Moorkamp  et  al.,  2007,  Jones  et  al.,  2009),  the  two 
studies  will  provide  some  independent  confirmation  of  one  another.  Finally,  the  two 
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techniques  can  be  used  in  conjunction  to  emphasize  the  strengths  and  mitigate  the 
weaknesses  of  the  other. 

The  teleseismic  migration  method,  involving  backpropagation  of  a  scattered 
wavefield  to  identify  the  scattering  points  in  the  subsurface,  is  particularly  suited 
to  the  illumination  of  sharp  geological  boundaries  (Rondenay  et  al.,  2009). 
Conversely,  MT  is  highly  sensitive  to  the  presence  of  conductive  fluid  and  melt 
phases,  and  we  have  shown  that  these  can  be  even  better  constrained  if  we  have 
a  priori  understanding  of  some  of  the  resistive  features  in  the  profile.  Figure  7 
shows  the  primary  inversion  models  from  the  teleseismic  migration  study  from 
Chapter  2  as  well  as  the  MT  study  from  Chapter  3.  While  both  of  these  models 
have  been  discussed  in  detail  in  their  respective  chapters,  we  look  at  them  now 
in  a  new  light.  Specifically,  we  want  to  consider  the  migrated  model  as  a  priori 
information  that  can  be  incorporated  into  a  new  MT  inversion  in  order  to  better 
constrain  the  conductive  features  in  the  latter. 

From  the  migrated  image  (fig.  7a),  we  are  able  to  clearly  identify  the  top  of  the 
descending  slab,  from  a  depth  of  ~20  km  at  the  coastline  (left  margin)  to  depths 
of  ~80  km  or  more.  We  can  identify  regions  where  serpentinization  has  occurred 
due  to  fluid  release  from  the  eclogitization  of  the  hydrated  basalts  of  the  upper 
crust  (A)  (Rondenay  et  al.,  2001,  Bostock  et  al.,  2002),  the  cool  nose  of  the  mantle 
wedge  (B)  (van  Keken  et  al.,  2010),  and  a  fluid  or  melt  phase  near  the  top  of  the 
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slab  at  depth  (C)  consistent  with  the  release  of  fluids  from  chlorite  harzburgite  in 
the  upper  mantle  and  subsequent  conversion  of  metastable  gabbros  in  the  lower 
crust  to  eclogite  (Hacker  et  al.,  2003,  John  and  Schenk,  2003).  We  can  also  clearly 
make  out  the  continental  Moho  in  the  east. 
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Figure  7:  Primary  images  from  the  teleseismic  migration  (fig.  7a)  and 
magnetotelluric  (fig.  7b)  studies  from  the  previous  chapters. 


In  the  MT  image  (fig.  7b),  we  are  not  able  to  make  out  the  contrast  at  the  top  of 
the  slab,  but  based  on  our  discussion  in  section  4,  we  recognize  that  the 
smoothness  constraint  would  cause  a  resistive  slab  to  retreat  from  the  conductors 
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at  A  and  C.  These  conductors  serve  to  confirm  our  interpretation  of  the  migrated 
image,  as  does  the  resistive  feature  at  B.  Additionally,  the  MT  image  shows  quite 
clearly  the  connection  between  the  conductive  fluid  melt  at  C  and  an  intense 
conductor  (D)  that  is  consistent  with  a  magma  chamber  feeding  the  nearby  Mt. 
Rainier  (Stockstill  et  al.,  2002,  Hildreth,  2007,  Pommier  and  LeTrong,  2007).  The 
conductive  pathway  between  these  two  features  also  provides  an  explanation  for 
the  disruption  of  the  continental  Moho  evident  in  the  migration  image. 

While  the  two  models  are  clearly  consistent  with  one  another,  the  objective  is 
to  improve  the  resolution  of  the  conductors  in  the  MT  model  using  the 
techniques  from  section  4.  Using  the  migration  model  to  define  the  upper 
boundary  of  the  descending  slab,  we  impose  a  resistor  into  the  starting  model 
and  a  tear  zone  over  the  top  few  km  of  the  imposed  resistor.  Additional 
inversions  were  run  with  a  tear  zone  at  the  continental  Moho  as  defined  by  the 
migration  model,  and  over  the  conductor  below  Mt.  Rainier,  but  these  did  not 
improve  the  images  and  were  not  incorporated  into  the  final  augmented  model 
(fig.  8). 

It  is  clear  that  the  constraints  on  the  conductors  are  somewhat  different  in  the 
augmented  model  as  compared  to  the  original  halfspace  inversion,  but  we  must 
recognize  that  we  have  traded  one  assumption  for  another  and  still  bear  the 
burden  of  demonstrating  that  the  new  model  is  more  accurate.  The  traditional 
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way  to  accomplish  this  is  to  look  at  the  rms  misfit  values  associated  with  each 
model.  The  rms  misfit  value  for  the  original  halfspace  model  was  3.08,  and  for 


Figure  8:  The  top  panel  shows  the  augmented  MT  inversion  model 
incorporating  a  priori  information  using  the  methods  developed  in  section  4.  In 
the  lower  panel,  a  thermal  profile  for  the  Cascades  (van  Keken,  2011)  is 
superimposed  on  the  model.  Note  that  the  thermal  model  assumes  a  shallower 
dip  angle  for  the  slab  in  the  eastern  part  of  the  model. 
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the  augmented  model  1.89  which  represents  a  significant  improvement.  The 
picture  becomes  even  clearer  when  we  look  at  a  plot  of  rms  values  at  different 
stations  (fig.  9). 


Figure  9:  Plot  of  rms  misfit  at  the  60  long  period  stations  for  the  TM/TE/tipper 
models  of  the  halfspace  inversion  from  Chapter  3,  and  the  augmented  inversion 
incorporating  the  a  priori  information  from  the  teleseismic  migration  using  the 
methods  developed  throughout  section  4  of  this  chapter.  The  overall  rms  values 
are  3.08  for  the  halfspace  inversion,  and  1.89  for  the  augmented  inversion. 


While  the  augmented  model  shows  improvement  in  fit  at  almost  every  station, 
the  effect  is  much  more  pronounced  between  stations  30  and  55,  which 
correspond  to  horizontal  distances  of  150-275  km  from  the  coast.  By  far  the  most 
striking  differences  between  the  two  models  between  these  markers  is  that  in  the 
augmented  model  the  conductor  is  more  intense  and  not  separated  from  the 
resistive  slab  by  any  significant  distance.  Taken  together,  this  strongly  suggests 
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that  the  augmented  model  is  more  accurate,  and  that  the  imposed  smoothness 
was  largely  responsible  for  the  higher  misfit  in  the  halfspace  model. 

A  second  way  to  evaluate  the  accuracy  of  the  models  is  by  generating 
synthetic  data  from  them,  and  inverting  the  synthetic  data  using  the  same 
parameters  that  were  used  for  the  original  inversions.  If  a  model  is  an  accurate 
representation  of  the  actual  subsurface,  then  the  synthetic  data  generated  from 
the  model  should  be  a  reasonable  copy  of  the  actual  data,  and  inversions  using 
the  synthetic  data  should  reasonably  reproduce  the  results.  Figure  10  shows  the 
results  of  using  the  original  halfspace  model  (fig.  10a)  as  the  forward  model.  The 
synthetic  data  is  then  used  to  generate  two  new  inversion  models,  one  using  the 
parameters  from  the  original  halfspace  model,  and  the  other  using  the 
parameters  from  the  augmented  model.  These  results  are  shown  in  fig.  10b  and 
fig.  10c  respectively. 

The  model  generated  from  the  halfspace  reproduces  many  of  the  features  of 
the  original  fairly  well,  but  there  are  some  significant  differences  as  regards  the 
conductor  rising  from  depth.  In  the  original  model,  the  resistivity  of  the 
conductor  rising  from  depths  of  80-100  km  is  about  7  Qm,  a  value  that  is 
maintained  from  depths  of  35-100  km.  For  the  reproduction,  this  resistivity  is 
only  found  to  depths  of  ~45  km,  and  at  100  km  has  risen  to  more  than  20  Qm. 
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Figure  10:  A  series  of  models  designed  to  test  the  accuracy  of  the  halfspace 
model.  Fig.  10a  is  used  as  the  starting  model  to  generate  the  synthetic  data.  Fig. 
10b  is  the  new  inversion  from  a  halfspace,  and  fig.  10c  is  the  new  inversion  using 
the  tools  from  the  augmented  inversion. 


The  dipping  resistive  feature  has  also  retreated  even  further,  suggesting  a 
significant  role  for  the  imposed  smoothness. 

The  augmented  reproduction  has  even  more  difficulty.  The  resistivity  at  the 
top  of  the  slab  is  ~20  Cm  compared  to  2-3  Qm  in  the  original  augmented  model. 
A  conductive  (3  Qm)  artifact  has  appeared  at  a  depth  of  35  km  within  the  column 
extending  from  the  slab.  If  the  original  halfspace  model  had  been  accurate,  then 
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imposition  of  the  dipping  resistor  and  tear  zone  would  have  concentrated  the 
conductivity  at  the  top  of  the  slab. 

Figure  11  shows  the  results  of  using  the  augmented  model  (fig.  11a)  as  the 
forward  model  for  the  synthetic  data.  The  new  reproduction  models  are  then 
generated  in  the  same  way  as  before,  with  the  first  using  parameters  from  the 
original  halfspace  model  and  the  second  from  the  augmented  model.  These 
results  are  shown  in  fig.  lib  and  fig.  11c  respectively. 

The  reproduction  in  the  halfspace  model  is  quite  good.  Since  we  know  what 
the  forward  model  looked  like,  we  can  see  the  effect  that  the  imposed 
smoothness  has  had  on  the  deep  conductor,  and  its  appearance  and  resistivity 
values  (7-8  Qm)  are  very  similar  to  what  we  see  in  the  original  halfspace  model. 
The  appearance  of  the  resistive  slab  is  also  very  similar,  again  suggesting  that  the 
imposed  smoothness  played  a  very  similar  role  in  this  inversion  as  it  did  for  the 
original  halfspace.  The  other  features  of  the  model  are  also  very  similar,  and  we 
conclude  that  for  this  inversion,  the  reproduction  is  quite  good. 

Turning  now  to  the  augmented  reproduction  we  see  that  the  most  conductive 
region  of  the  column  is  at  the  top  of  the  slab,  with  resistivities  approaching  5  Qm. 
While  there  is  some  concentration  of  conductivity  at  35  km,  the  overall  shape  of 
the  conductor  is  consistent  with  the  original  image.  The  other  features  in  the 
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Figure  11:  A  series  of  models  designed  to  test  the  accuracy  of  the  augmented 
model.  Fig.  11a  is  used  as  the  starting  model  to  generate  the  synthetic  data.  Fig. 
lib  is  the  new  inversion  from  a  halfspace,  and  fig.  11c  is  the  new  inversion  using 
the  tools  from  the  augmented  inversion. 


model  are  also  quite  accurately  reproduced.  While  the  reproduction  is  certainly 
imperfect,  the  improvement  relative  to  the  models  generated  from  the  halfspace 
starting  model  provides  a  stark  contrast,  and  we  can  only  conclude  that  the 
augmented  model  is  a  superior  representation  of  the  conductivity  structure  of 
the  subsurface. 
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6.  Implications 


Figure  12  shows  the  interpreted  seismic  migration  image  from  the  study 
conducted  in  Chapter  2,  along  with  the  augmented  MT  model  from  this  chapter. 
The  interpretation  of  both  the  seismic  image  and  the  original  MT  image  was 
conducted  in  detail  in  the  earlier  chapters,  so  we  will  not  repeat  all  of  the 
previous  analysis.  The  augmented  MT  image  does  however  provide  better 
constraints  for  some  of  the  conductive  features,  and  the  implications  pertaining 
to  them  needs  to  be  addressed. 

The  conductor  (A)  that  appears  to  be  in  contact  with  the  slab  at  depths  of  40 
km  which  had  a  resistivity  of  80  Qm  in  the  original  value  is  now  ~  40  Qm,  more 
in  line  with  what  has  been  recorded  in  other  subduction  settings  (Brasse  et  al., 
2009,  Evans  et  al.,  submitted),  where  it  has  been  interpreted  as  migrating  fluids 
that  have  been  released  from  the  subducting  oceanic  crust.  This  resistivity  is 
consistent  with  a  0.5%wt  of  aqueous  fluid  given  good  connectivity;  a  tradeoff 
between  higher  fluid  content  and  lower  connectivity  also  produces  an  acceptable 
result  (Pommier  and  LeTrong,  2007). 

This  conductor  can  also  be  seen  to  extend  up-dip  along  the  top  of  the  slab, 
supporting  an  earlier  interpretation  of  seismic  results  indicating  the  presence  of 
pressurized  channels  related  to  episodic  tremor  and  slip  (ETS)  (Abers  et  al.,  2009, 
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Audet  et  al.,  2010).  This  feature  does  not  appear  to  be  present  along  the  EMSLAB 


transect  (Evans  et  al.,  submitted),  where  ETS  is  also  absent. 


50  100  150  200  250 

Horizontal  distance  (km) 


Figure  12:  Final  seismic  migration  and  MT  augmented  inversion  models  for  the 
CAFE  experiment. 


The  deeper  conductor  (B)  has  low  resistivities  (~  3  Qm)  along  the  top  of  the 
slab,  but  sensitivity  tests  on  this  feature  suggest  that  the  resistivity  of  this  feature 
could  be  as  high  as  8  Qm  (see  Appendix).  In  practice,  at  these  depths  even  this 
higher  value  is  not  so  easy  to  achieve.  Based  on  an  adaptation  of  the  Vogel- 
Fulcher-Tammann  equation  (eq.  1)  (Ni  et  al.,  2011)  and  the  HS+  bounds  (eq.  2) 
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(Hashin  and  Shtrikman,  1962),  we  can  calculate  the  bulk  conductivity  based  on 
only  temperature,  weight  percent  water,  and  volume  fraction  of  the  melt. 


logo-  =  2.172  — 


860.82-206.46  W^2 
T- 1146.8 


(1) 


W  is  the  H2O  content  in  wt%,  T  is  the  temperature  in  K,  and  a  is  the  electrical 
conductivity  in  S/m. 


O  =  o2 


3  02 — F2  (er2 — 04). 


(2) 


F2  is  the  melt  fraction,  02  is  the  conductivity  in  S/m  of  the  more  conductive  phase, 
and  a  is  the  bulk  conductivity. 


Table  1  shows  the  values  of  bulk  resistivities  calculated  from  equations  1  and 
2  given  a  temperature  of  1200  C,  a  background  conductivity  of  .01  S/m,  and 
varying  wt%  H2O  from  0-5%,  and  melt  fraction  from  1-4%. 

The  conductor  appears  to  be  in  contact  with  the  slab  from  depths  of  75-105 
km,  corresponding  to  the  range  of  dehydration  for  chlorite  harzburgite  in  the 
mantle  (Hacker  et  al.,  2003)  and  the  transition  of  metastable  gabbro  to  eclogite 
(John  and  Schenk,  2003)  for  Cascadia  (van  Keken  et  al.,  2011).  This  conductor 
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Weight  percent  FRO 

0  wt.% 

1  wt.% 

2  wt.% 

3  wt% 

4  wt% 

Melt  fraction 

1% 

81 

51 

36 

26 

20 

2% 

69 

34 

22 

15 

11 

3% 

60 

26 

16 

11 

7 

4% 

52 

20 

12 

8 

6 

5% 

47 

17 

10 

7 

5 

Table  1:  Bulk  resistivity  (Cm)  of  fluid/melt  phase  at  the  top  of 


the  slab  for  a  selection  of  melt  fraction  and  water  content 


extends  upward  with  a  resistivity  of  ~  7-10  Cm  until  it  connects  with  the  shallow 
conductor  (C)  of  resistivities  between  3-4  Cm  that  we  interpret  to  be  a  magma 
chamber  that  feeds  Mt.  Rainier. 

While  the  smoothness  imposed  on  standard  MT  inversions  is  a  desirable 
feature  that  stabilizes  the  inversion  and  minimizes  the  introduction  of  structure 
that  is  not  required  by  the  data,  the  penalty  it  applies  to  abrupt  transition  can 
have  a  detrimental  effect  on  resolving  certain  features.  In  this  chapter,  we 
examine  how  this  assumption  of  smoothness  can  impact  the  final  inversion 
model,  methods  such  as  the  imposition  of  resistivities  and  tear  zones  that  can  be 
used  to  mitigate  those  effects,  and  ways  to  evaluate  which  set  of  assumptions 
provides  a  more  accurate  inversion  model.  Finally,  we  applied  the  lessons  to  the 
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CAFE  dataset,  resulting  in  a  more  accurate  model,  particularly  with  respect  to 
the  constraints  for  the  conductive  features  in  the  model. 
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Appendices 


The  appendices  are  presented  as  a  series  of  figures  designed  to  provide  a 
deeper  look  into  the  mechanics  of  the  work  done  while  constructing  this  thesis. 
The  first  appendix  supports  the  seismic  work,  presenting  the  steps  involved  in 
the  preprocessing  of  the  seismic  data  and  the  migrated  images  generated  for  each 
of  the  64  events  that  were  used  in  the  final  images.  The  second  appendix 
supports  the  MT  work,  presenting  pseudo-sections  derived  from  the  original 
(processed)  data,  difference  pseudo-sections  comparing  the  augmented  model 
with  the  standard  model,  and  an  inversion  model  generated  in  the  absence  of  a 
smoothing  parameter  to  support  the  theory  presented  in  Chapter  4. 

Additionally,  some  of  the  most  important  sensitivity  tests  are  included,  such  as 
the  test  for  robustness  of  the  vertical  conductive  corridor  that  connects  the  two 
most  prominent  conductors  and  the  allowable  range  test  for  the  conductivities  at 
the  top  of  the  slab.  Finally,  some  synthetic  models  are  included  that  explore  the 
ramifications  of  imposing  features  and/or  tear  zones  based  on  faulty  a  priori 
information. 


153 


5.1  Figures  supporting  seismic  migration  work 
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Figure  1:  After  rotation  of  the  wavefield  and  application  of  the  free  surface  transfer 
matrix,  the  next  step  in  the  pre-processing  involved  the  manual  inspection  of  the 
waveform  for  each  of  the  41  stations  for  approximately  250  prospective  events.  The  first 
three  panels  of  this  figure  show  typical  waveforms,  taken  from  stations  N040,  N015,  and 
S010  for  event  20082541308  (the  event  designator  consists  of  the  year,  Julian  Date,  hour, 
and  minute  of  the  event).  A  prominent  peak  was  then  selected  that  could  be  readily 
identified  for  each  station  (shown  in  the  figure  as  the  vertical  red  line).  Manual 
inspection  of  the  waveforms  was  also  conducted  in  the  frequency  domain  (the  fourth 
panel  shows  station  S015  for  the  same  event).  In  general,  the  frequencies  were  filtered  to 
pass  those  between  0.03  and  0.5  Hz,  but  some  stations,  particularly  those  with  excessive 
low  frequency  noise,  required  additional  filtering. 
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Figure  2:  The  initial  picks  of  prominent  peaks  shown  in  figure  1  ensure  that  the  start  times  for 
each  event  are  reasonably  well  aligned  between  stations.  This  alignment  is  then  refined  using  a 
multi-channel  cross-correlation  technique-  figure  2a  shows  the  aligned  waveforms  for  39 
stations  for  event  20082541308.  The  next  step  in  the  pre-processing  is  the  application  of 
eigenimage  decomposition  to  identify  the  correlated  parts  of  these  waveforms,  under  the 
assumption  that  these  correlated  parts  provide  a  good  estimate  of  the  source  waveform.  The 
source  waveform  for  event  20082541308  is  shown  in  figure  2b.  The  next  step  is  to  deconvolve 
the  estimated  incident  wavefield  from  the  components  of  the  scattered  wavefield  to  obtain  the 
normalized  scattered  wavefield  (corresponding  to  a  source  impulse).  In  order  to  make  this 
deconvolution  stable,  a  "waterlevel"  proportional  to  noise  is  used  in  the  denominator  of  the 
deconvolution.  For  this  experiment,  we  used  a  station  specific  waterlevel  (developed  by  Fred 
Pearce  at  MIT)  so  that  noisy  stations  would  not  cause  loss  of  data  for  all  of  the  stations  for  a 
given  event.  Figure  2c  shows  the  waterlevel  plot  for  event  20082541308. 
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Figure  3:  The  deconvolution  of  the  estimated  source  wavefield  results  in  a  three-component 
estimate  of  the  scattered  wavefield  for  a  specific  event.  These  scattered  wavefields  are  then 
migrated  back  to  find  the  scattering  points  in  the  subsurface.  This  figure  shows  the  scattered 
wavefields  for  event  20082541308. 
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Figure  4:  The  individual  composite  images  for  each  of  the  64  events  that  were  used  to  generate 
the  primary  composite  image.  The  modes  used  for  each  of  these  images  are  given  in  the 
supplemental  table  to  Chapter  2. 


200  100 
distance  along  projection  line  (km) 


20081970326 


20082010239 


20082521852 


20082541308 


20082400135 


20082392100 


20082051526 


20082010927 


164 


5.2  Additional  figures  supporting  magnetotellurics  work 


u 

01 


T3 

O 

'u 

01 

Q. 


O 

O 

o 

o 

o 


Apparent  Resistivity  (ohm  m) 


o 

o 

o 

o 


o 

o  o 

o  o 


o 


Distance  from  coast  (km) 

Figure  5:  TM  mode  pseudo  section  for  the  CAFE  MT  data.  The  upper  figure  shows  apparent 
resistivity  and  the  lower  figure  shows  phase.  All  of  the  pseudo  section  models  (figures  5-8)  are 
limited  horizontally  to  correspond  with  the  surface  covered  by  the  CAFE  MT  array. 
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Figure  6:  TE  mode  pseudo  section  for  the  CAFE  MT  data.  The  upper  figure  shows  apparent 
resistivity  and  the  lower  figure  shows  phase. 
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Figure  7:  Pseudo  sections  can  also  be  generated  from  model  derived  synthetic  data,  and  models 
can  be  compared  using  difference  pseudo  sections.  This  figure  shows  the  difference  pseudo 
section  for  the  TM  mode  between  the  augmented  and  standard  models  from  Chapter  3.  Both 
models  fit  the  data  very  well  in  the  TM,  and  the  differences  between  them  are  quite  small. 


167 


Apparent  Resistivity  (%  difference) 

8  °  ° 

(j  oo 

M  I  I  I  I  I  1  I  I  I  II  I 


Distance  from  coast  (km) 

Figure  8:  This  figure  shows  the  difference  pseudo  section  for  the  TE  mode  between  the 
augmented  and  standard  models  from  Chapter  3.  The  differences  here  are  significantly  more 
pronounced  that  they  were  for  the  TM  mode,  which  is  not  unexpected  given  the  poor  TE  fit  for 
the  standard  model.  It  is  clear  again  here  that  the  most  important  differences  occur  between 
100-1000  seconds  roughly  200  km  from  the  coast,  consistent  with  the  fluid  release  from  the  slab 
at  depth. 
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_ Horizontal  distance  (km) _ 

Figure  9:  Figure  9a  shows  an  example  of  the  sort  of  model  that  is  ruled  out  by  the  use  of  a 
smoothness  driven  stabilizing  functional  in  an  MT  inversion.  Figure  9b  shows  an  inversion 
model  generated  from  the  CAFE  MT  data  for  which  the  stabilizing  parameter  has  been  collapsed 
by  imposing  a  tear  zone  over  the  entire  cross-section.  The  tau  parameter  was  also  set  to  the 
minimum  0.01  so  as  to  have  the  smallest  possible  impact  on  the  direction  of  the  inversion. 
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Figure  10:  Sensitivity  analysis  is  generally  conducted  on  all  of  the  primary  features  of  the 
inversion  models,  as  discussed  in  Chapters  3  and  4.  This  figure  shows  a  typical  example  in  which 
the  robustness  of  the  connecting  corridor  between  the  lower  conductor  associated  with  the  slab 
and  the  upper  conductor  associated  with  Mt.  Rainier.  The  primary  augmented  model  from 
Chapter  4  is  shown  in  figure  10a.  In  figure  10b,  the  vertical  conductive  corridor  has  been 
removed  and  replaced  by  a  region  with  a  conductivity  of  100  Dm.  The  inversion  was  restarted 
from  this  model,  resulting  in  the  model  shown  in  figure  10c.  It  is  clear  that  the  conductive 
corridor  is  being  regenerated  in  this  model,  suggesting  that  it  is  a  robust  feature. 
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Figure  11:  Another  type  of  sensitivity  analysis  consists  of  varying  the  resistivity  of  a  feature  and 
evaluating  the  effect  of  the  variations  on  misfit.  This  figure  shows  one  of  a  series  of  tests 
conducted  on  the  deep  conductor  in  the  augmented  model  to  constrain  its  allowed  conductivity. 
Figure  11a  shows  the  primary  augmented  model  from  Chapter  4.  Figure  lib  shows  the  new 
starting  model  for  which  the  minimum  resistivity  of  the  conductor  is  10  Dm.  Figure  11c  shows 
the  inversion  result:  the  inversion  was  able  to  better  the  original  rms  misfit  value  while  the 
minimum  resistivity  of  the  conductor  dropped  to  8  Dm.  This  series  of  sensitivity  tests  constrains 
the  lowest  resistivity  of  the  conductor  rising  from  the  slab  to  be  between  2  and  8  Dm,  with 
implications  for  its  possible  composition. 
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Figure  12:  Part  of  the  task  of  evaluating  synthetic  models  is  to  investigate  how  incorrect 
assumptions  can  influence  the  results.  Figure  12a  shows  the  synthetic  forward  model  from 
figure  2a  in  Chapter  4,  where  a  dipping  resistive  (9000  Dm)  feature  is  in  close  proximity  to  two 
relatively  thin  conductors  (3  Qm).  This  model  was  used  to  generate  a  set  of  synthetic  data. 
Figure  12b  shows  a  starting  model  for  an  inversion  series  for  which  particularly  poor 
assumptions  were  made  regarding  the  location  of  the  resistive  feature.  Figure  12c  shows  the 
inversion  model  generated  from  the  poorly  selected  starting  model.  It  is  clear  that  the  initial 
assumptions  have  been  largely  corrected  by  the  inversion,  an  impression  that  is  confirmed  by 
comparing  this  model  to  the  halfspace  inversion  model  shown  in  figure  2c  in  Chapter  4. 
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Figure  13:  These  models  also  use  the  synthetic  data  generated  from  the  forward  model  in  figure 
12a,  as  well  as  the  poor  assumptions,  but  this  time  using  tear  zones.  Figure  13a  shows  the 
starting  model  as  a  100  Qm  half  space  with  tear  zones  (but  no  resistive  features).  The  inversion 
model  associated  with  this  starting  model  is  shown  in  figure  13b.  While  the  locations  of  the 
conductors  are  fairly  accurately  resolved,  the  interface  between  the  shallow  conductor  and  the 
resistive  feature  incorrectly  follows  the  imposed  tear  zone.  Figure  13c  shows  a  starting  model 
for  which  both  the  resistive  features  and  their  corresponding  tear  zones  are  included,  and  the 
associated  inversion  model  is  shown  in  figure  13d.  In  this  case,  the  inversion  model  is 
reasonably  accurate  away  from  the  tear  zone,  but  retains  the  incorrect  assumption  implicit  in 
the  starting  model  within  the  tear  zone.  The  eastern  resistive  feature  is  much  less  prominent  in 
the  inverted  model,  but  still  present  as  well.  Figures  12  and  13  demonstrate  the  importance  of 
using  resistive  features  first  without  tears  in  order  to  evaluate  the  accuracy  of  the  assumptions 
about  location;  this  was  a  recurring  theme  in  our  analyses. 
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Figure  14:  Figure  14  shows  the  impact  of  an  assumption  in  which  the  top  of  the  resistive  clipping 
feature  is  placed  approximately  10  km  too  high  in  the  starting  model  (again  based  on  synthetic 
data  generated  from  the  forward  model  in  figure  12a).  Figure  14a  is  the  starting  model  showing 
the  assumed  location  of  the  resistive  feature.  The  associated  inversion  model  in  figure  14b  does 
a  nice  job  relocating  the  dipping  resistor  to  its  proper  location,  though  for  reasons  explained  in 
Chapter  4  does  not  constrain  the  conductors  very  well.  Figure  14c  shows  a  starting  model  with 
the  same  imposed  resistive  feature,  but  with  the  additional  imposition  of  a  tear  zone  along  the 
top  of  the  resistor.  This  causes  the  inversion  (figure  14c)  to  incorrectly  hold  the  resistor  in  the 
position  of  the  tear. 
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Figure  15:  Figure  14  shows  the  impact  of  an  assumption  in  which  the  top  of  the  resistive  clipping 


feature  is  placed  approximately  10  km  too  low  in  the  starting  model  (again  based  on  synthetic 
data  generated  from  the  forward  model  in  figure  12a).  The  inversion  in  figure  14b  moves  the 
resistive  feature  shallower,  though  it  also  exhibits  "retreat"  from  the  conductive  feature  as 
described  in  Chapter  4.  The  starting  model  in  figure  14c  includes  an  imposed  tear  zone,  and  the 
associated  inverted  model  is  shown  in  figure  14d.  In  this  case,  the  resistive  feature  passes  the 
tear  zone  boundary  and  extends  upwards.  These  kinds  of  analyses  show  that  examining  starting 
models  with  combinations  of  imposed  features  and  tear  zones,  we  can  determine  with  some 
confidence  both  the  accuracy  of  our  assumptions  and  the  manner  in  which  they  might  not  be 
correct. 
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